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Abstract 

The  quantification  of  uncertainties  is  critical  when  systems  are  nonlinear  and  have  uncer¬ 
tain  terms  in  their  governing  equations  or  are  constrained  by  limited  knowledge  of  initial 
and  boundary  conditions.  Such  situations  are  common  in  multiscale,  intermittent  and  non- 
homogeneous  fluid  and  ocean  flows.  The  Dynamically  Orthogonal  (DO)  field  equations 
provide  an  efficient  time-dependent  adaptive  methodology  to  predict  the  probability  den¬ 
sity  functions  of  such  flows.  The  present  work  derives  efficient  computational  schemes  for 
the  DO  methodology  applied  to  unsteady  stochastic  Navier-Stokes  and  Boussinesq  equa¬ 
tions,  and  illustrates  and  studies  the  numerical  aspects  of  these  schemes.  Semi-implicit 
projection  methods  are  developed  for  the  mean  and  for  the  orthonormal  modes  that  define 
a  basis  for  the  evolving  DO  subspace,  and  time-marching  schemes  of  first  to  fourth  order 
are  used  for  the  stochastic  coefficients.  Conservative  second-order  finite-volumes  are  em¬ 
ployed  in  physical  space  with  Total  Variation  Diminishing  schemes  for  the  advection  terms. 
Other  results  specific  to  the  DO  equations  include:  (i)  the  definition  of  pseudo-stochastic 
pressures  to  obtain  a  number  of  pressure  equations  that  is  linear  in  the  subspace  size  in¬ 
stead  of  quadratic;  (ii)  symmetric  Total  Variation  Diminishing-based  advection  schemes 
for  the  stochastic  velocities;  (iii)  the  use  of  generalized  inversion  to  deal  with  singular 
subspace  covariances  or  deterministic  modes;  and  (iv)  schemes  to  maintain  orthonormal 
modes  at  the  numerical  level.  To  verify  the  correctness  of  our  implementation  and  study 
the  properties  of  our  schemes  and  their  variations,  a  set  of  stochastic  flow  benchmarks  are 
defined  including  asymmetric  Dirac  and  symmetric  lock-exchange  flows,  lid-driven  cav¬ 
ity  flows,  and  flows  past  objects  in  a  confined  channel.  Different  Reynolds  number  and 
Grashof  number  regimes  are  employed  to  illustrate  robustness.  Optimal  convergence  un¬ 
der  both  time  and  space  refinements  is  shown  as  well  as  the  convergence  of  the  probability 
density  functions  with  the  number  of  stochastic  realizations. 
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1.  Introduction 

Quantifying  uncertainty  is  becoming  increasingly  important  in  many  sci¬ 
entific  and  engineering  applications.  This  is  in  part  because  the  accuracy  of 
an  answer  is  now  often  as  critical  as  the  answer  itself.  As  a  result,  stochastic 
modeling  methods  for  uncertainty  quantification  are  being  developed.  Our 
present  motivation  is  uncertainty  quantification  for  computational  fluid  dy¬ 
namics  (CFD)  applications,  specifically  in  the  context  of  providing  realistic 
predictions  of  ocean  fields.  In  ocean  dynamics,  it  is  challenging  to  model 
multi-scale,  intermittent,  non-stationary  and  non- homogeneous  uncertain¬ 
ties.  Already  a  single  evaluation  of  an  ocean  model  is  costly  and  straight¬ 
forward  stochastic  modeling  methods  are  prohibitively  expensive  [38,  48], 
particularly  when  dealing  with  longer-term  unsteady  nonlinear  dynamics. 
Fortunately,  the  recently  developed  Dynamically  Orthogonal  (DO)  field  equa¬ 
tions  [52,  51,  53]  provide  an  efficient,  tractable  means  for  uncertainty  pre¬ 
diction  in  large-scale  CFD  and  ocean  applications.  While  these  DO  equa¬ 
tions  have  been  solved  numerically  using  a  simple  finite-difference  scheme, 
the  specific  properties  of  the  DO  equations  warrant  novel  integration  and 
discretization  schemes.  Hence,  our  present  goals  are  to  derive  efficient  com¬ 
putational  schemes  for  the  DO  methodology  applied  to  unsteady  stochastic 
Navier-Stokes  and  Boussinesq  equations,  and  to  illustrate  and  study  the  nu¬ 
merical  aspects  of  these  schemes. 

Stochastic  modeling  approaches  can  be  categorized  as  either  non- intrusive 
or  intrusive.  Non-intrusive  approaches  have  the  advantage  that  the  deter¬ 
ministic  version  of  a  model  can  be  used  to  generate  an  ensemble  of  sample 
solutions  from  which  the  statistics  can  be  calculated.  The  disadvantage  is 
that  a  large  number  of  samples  are  often  required,  leading  to  large  computa¬ 
tional  costs.  Intrusive  approaches  require  the  update  of  existing  codes  or  the 
development  of  a  new  code,  where  the  resulting  system  of  equations  can  be 
larger  than  the  original,  deterministic  system.  However,  intrusive  methods 
are  usually  computationally  less  expensive  than  non-intrusive  methods,  and 
the  statistics  are  explicitly  available.  Below  we  give  a  brief  review  of  existing 
methods,  focusing  on  their  computational  aspects.  For  more  complete  refer¬ 
ences  and  reviews  we  refer  for  e.g.  to:  Ghanem  and  Spanos  [18],  Kloeden  and 
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Platen  [26],  Doucet  et  al  [9],  Mathelin  et  al  [45],  Eldred  et  al  [10],  Jakeman 
and  Roberts  [24],  Najm  [49],  Xiu  [67,  68],  Le  Matre  and  Knio  [30]. 

The  non-intrusive  Monte-Carlo  method  provides  access  to  the  full  statis¬ 
tics  of  the  problem.  Its  computational  cost  does  not  strictly  depend  on  the 
size  of  the  system,  but  more  on  the  number  of  truly  independent  random 
variables,  and  convergence  rates  are  often  proportional  to  the  square  root 
of  the  number  of  samples.  The  efficiency  can  be  improved  for  example  by 
using  more  elaborate  Monte-Carlo  schemes  [e.g.  9],  including  particle  filters 
or  mixtures  of  weighted  kernels,  e.g.  Gaussian  kernels  [5,  47].  Nonetheless, 
a  large  number  of  function  evaluations  are  needed  due  to  the  slow  (square 
root)  convergence,  which  can  limit  accuracy  in  large-scale  applications. 

The  Polynomial  chaos  expansion  (PCE),  pioneered  by  Ghanem  and  Spanos 
[18]  and  based  on  the  theory  by  Wiener  [65,  2,  66],  has  become  popular 
because  it  can  represent  and  propagate  large  uncertainties  through  com¬ 
plex  models.  Both  non-intrusive  [e.g.  64,  27,  23,  10]  and  intrusive  [e.g. 
7,  28,  62,  44,  49]  versions  have  been  employed,  but  both  can  suffer  from 
the  curse  of  dimensionality.  That  is,  the  PCE  of  a  dynamical  model  scales 
as  (p+f ,  where  p  is  the  largest  degree  of  the  polynomials  used,  and  s  is  the 
number  of  independent  random  variables.  For  problems  that  require  large 
p  (e.g.  non-Gaussian)  and  large  s  (e.g.  large  ocean  or  fluid  simulations), 
the  number  of  function  evaluations  (non-intrusive),  or  the  size  of  the  sys¬ 
tem  of  equations  (intrusive),  can  quickly  become  prohibitive.  For  s  larger 
than  p,  the  storage  of  a  PCE  scales  as  0(sp )  while  its  computational  cost 
for  Navier-Stokes  flows  scales  as  0(s2p )  due  to  the  quadratic  nonlinearity  of 
advection  terms.  The  large  cost  of  these  methods  have  prompted  the  use 
of  non-Gaussian  random  variables  [50],  the  development  of  generalized  PCE 
[69]  to  speed  up  convergence  in  the  polynomial  degree  (i.e.  reduce  p),  and 
the  development  of  adaptive  schemes  that  only  evaluate  the  necessary  terms 
in  the  PCE  (Li  and  Ghanem  [40]). 

PCEs  have  been  successful  in  many  CFD  applications.  In  the  case  of  un¬ 
steady  incompressible  fluid  dynamics,  Le  Maitre  et  al  [29]  used  a  PCE  scheme 
to  study  mixing  in  a  two-dimensional  (2D)  micro-channel  and  improved  the 
efficiency  of  their  solution  scheme  by  decoupling  the  velocity-pressure  equa¬ 
tions  using  a  projection  method.  Wan  and  Karniadakis  [63]  studied  the  long 
term  behavior  of  a  generalized  PCE  first  using  the  ID  advection  equation 
and  then  a  2D  noisy  flow  past  a  circular  cylinder.  The  authors  showed  that 
multi-element  generalized  PCEs  can  significantly  improve  the  accuracy  for 
long  time  integration,  but  caution  that  the  cost  for  large  s  remain  high.  Other 
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applications  include  fluid-structure  interactions  [e.g.  70],  turbulence  [e.g.  43], 
aerodynamics  [e.g.  56],  and  ocean  acoustics  [e.g.  14].  Other  examples  are 
also  studied  in  the  above  references. 

Motivated  by  the  multi-scale,  intermittent  and  non-homogeneous  uncer¬ 
tain  ocean  fields,  the  Error  Subspace  Statistical  Estimation  (ESSE)  method 
was  developed  [37,  31,  32],  It  uses  a  Karhunen-Loeve  (KL)  expansion  [25,  41] 
but  with  time  varying  and  adaptive  basis  functions.  This  generalized  KL  ex¬ 
pansion  is  initialized  by  a  multi-scale  scheme  [33]  and  evolved  using  stochas¬ 
tic,  data-assimilative  and  adaptive,  Monte-Carlo  ensemble  schemes.  The 
computational  cost  of  ESSE  predictions  scales  as  the  size  q  of  the  Monte- 
Carlo  ensemble  required,  i.e.  0(q).  This  ensemble  size  q  varies  with  time 
and  is  linked  to,  but  a  bit  larger  than,  the  evolving  size  of  the  error  subspace 
itself  s  (which  gives  the  storage  cost  scaling  in  0(s)). 

The  DO  equations  for  dynamically  evolving  stochastic  fields  [52,  51]  were 
derived  specifically  to  optimally  approximate  the  Fokker-Planck  equations 
and  capture  the  dominant  stochastic  subspace  while  being  computationally 
tractable.  The  DO  methodology  also  starts  from  a  truncated  generalized 
Karhunen-Loeve  expansion  but  derives  the  governing  equations  for  the  mean, 
the  modes  and  their  coefficients.  In  this  derivation,  a  key  condition  is  im¬ 
posed:  the  rate-of-change  of  the  stochastic  subspace  is  dynamically  orthog¬ 
onal  to  the  subspace  itself.  The  DO  subspace  basis,  i.e.  the  DO  modes,  as 
well  as  the  probability  density  functions  (pdfs),  i.e.  the  stochastic  coefficients, 
thus  evolves  only  according  to  the  dynamics  of  the  system.  This  renders  the 
DO  decomposition  efficient  and  limits  divergence  issues.  Consequently,  the 
computational  scaling  is  only  dependent  on  the  number  of  random  variables, 
s,  even  when  dealing  with  non-Gaussian  processes:  specifically,  the  storage 
scales  as  O(s)  and  computational  cost  as  0(s2)  for  Navier-Stokes  equations. 
The  size  s  is  in  general  a  function  of  time  so  as  to  adapt  to  the  dynamically 
evolving  uncertainties  and  boundary  conditions  [53].  The  DO  methodology, 
first  solved  numerically  in  Sapsis  and  Lermusiaux  [52],  has  been  applied  to 
several  Navier-Stokes  flows  and  their  stochastic  dynamics  has  been  stud¬ 
ied,  including  mean-mode  and  mode-mode  energy  transfers  for  2D  flows  and 
heat  transfers  [51,  54],  However,  the  DO  method  has  not  yet  been  applied 
to  Boussinesq  flows  relevant  for  ocean  studies,  and  the  numerical  challenges 
of  the  DO  decomposition  for  the  stochastic  Navier-Stokes  and  Boussinesq 
equations  have  not  be  thoroughly  examined  nor  resolved.  This  explains  the 
need  for  the  present  study. 

In  what  follows,  the  equations  for  incompressible  stochastic  Navier-Stokes 
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and  Boussinesq  equations  are  given,  and  their  DO  decomposition  is  outlined 
(§2).  In  §3,  the  discretization  in  time  is  developed,  discussing  explicit  and  im¬ 
plicit  schemes.  Importantly,  we  first  combine  DO  pressure  terms  and  define 
new  “pseudo-stochastic  pressures”  in  the  original  continuous  DO  equations. 
With  this  new  definition,  given  that  the  pressure  Poisson  equations  and  the 
corresponding  matrix  inversions  often  dominate  the  computational  cost,  the 
cost  of  DO  integration  schemes  is  substantially  reduced:  instead  of  scaling 
as  0(s2)  [53],  it  scales  as  O(s)  (as  long  as  the  solution  of  the  pressure  domi¬ 
nates  the  cost  of  the  scheme).  The  time  integration  schemes  are  then  derived. 
For  the  mean  and  the  modes,  we  employ  projection  methods  [19],  outlining 
schemes  of  first  and  second  order.  For  the  stochastic  coefficients,  we  obtain 
several  time-marching  schemes  of  first  to  fourth  order,  and  briefly  discuss 
their  extension  to  the  cases  of  additive  and  multiplicative  random  forcing. 
The  discretizations  of  the  physical  space  and  stochastic  subspace  are  given  in 
§4.  For  the  former,  the  discretization  of  diffusion  operators  is  straightforward. 
That  of  advection  operators  requires  special  attention:  since  the  determinis¬ 
tic  DO  modes  have  arbitrary  signs,  how  to  apply  upwinding  in  the  advection 
scheme  to  ensure  the  total  variation  diminishing  property  is  a  key  question  we 
investigate.  For  the  stochastic  subspace,  a  number  of  possible  discretizations 
are  outlined,  including  the  direct  Monte-Carlo  scheme.  In  §5,  the  questions 
of  how  to  deal  with  singular  covariances  and  how  to  maintain  orthonormal 
modes  in  presence  of  numerical  round-off  and  truncation  errors  are  discussed. 
For  the  applications  in  §6,  a  set  of  benchmarks  are  defined  and  utilized  to 
illustrate  the  properties  of  our  DO  numerics  scheme.  Specifically,  a  verifica¬ 
tion  benchmark  based  on  an  asymmetric  Dirac-stochastic  lock-exchange  flow 
is  defined  and  used  to  test  the  implementation.  A  symmetric  stochastic  lock- 
exchange  is  then  employed  to  evaluate  the  new  advection  schemes  for  DO 
modes.  The  spatial  and  temporal  convergence  is  studied  with  a  stochastic 
lid-driven  cavity  flow.  The  discretization  of  the  stochastic  coefficients  is  ex¬ 
amined  using  a  flow  over  a  square  cylinder  in  a  confined  channel.  Each  flow 
benchmark  is  purposely  chosen  to  be  different  in  part  to  illustrate  the  robust¬ 
ness  of  our  DO  numerics.  We  also  expect  that  such  benchmarks  can  serve  as 
standard  tests  for  future  schemes  and  implementations.  Lastly,  conclusions 
and  discussions  are  in  57. 
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2.  Stochastic  Discretization  of  Boussinesq  Equations 

This  section  defines  the  differential  equations  that  we  solve  and  then 
briefly  outlines  the  stochastic  DO  methodology.  The  temporal  and  spatial 
discretizations  are  derived  in  the  subsequent  section. 

The  deterministic  components  of  the  partial  differential  equations  (PDEs) 
that  we  solve  on  a  domain  T>  are  non-dimensional  Boussinesq  equations2,  in 
the  same  form  as  in  Hartel  et  al  [21], 

V  •  u  =  0,  x£D, 

W“7bv2u  =  ~v'(uu)-vp+pe!’’  xeV'  (i) 

%  ~  = -v ' (up)-  xeP- 

The  non-dimensional  variables  are:  u(x, t)  =  [u,v,w\,  the  velocity  in  3D; 
p(x,  t),  the  density;  and,  p(x,  t),  the  pressure.  The  vector  e9  is  a  unit-vector 
in  the  direction  of  gravity,  (x,  t)  are  the  non-dimensional  space  and  time 
variables,  Gr  =  9  ^  /  is  the  Grashof  number  which  is  the  ratio  of  buoyancy 
forces  to  viscous  forces,  Sc  =  u / K  is  the  Schmidt  number  which  is  the 
ratio  of  kinematic  viscosity  z>  to  molecular  diffusivity  K  for  the  density  field, 
g'  =  js  the  reduced  gravity,  and  h  is  the  vertical  length-scale.  In 

what  follows,  we  denote  the  total  dynamical  rate-of- change  in  the  prognostic 
eqns.  (1)  for  velocity  and  density  by  £u  and  £p,  respectively,  i.e.  =  £u 
and  | ^  =  £p. 

The  latter  prognostic  equation  for  density  originates  from  the  thermody¬ 
namic  energy  equation  and  an  equation  of  state  (it  arises  from  another  form 
of  the  Boussinesq  approximation  frequently  used  in  ocean  modeling  which 
retains  the  temperature  and  salinity  fields  as  state  variables,  e.g.  [6,  20]).  We 
emphasize  that  for  problems  without  density-driven  flows,  \jGr  =  Re ,  that 
is,  the  square  root  of  the  Grashof  number  is  the  Reynolds  number.  The  ap¬ 
proach  and  numerical  schemes  that  we  derive  in  this  manuscript  are  directly 
applicable  to  the  Navier-Stokes  equations. 

We  are  interested  in  solving  the  above  equations  in  their  stochastic  form, 
due  to  the  presence  of  uncertainties.  We  thus  introduce  the  set  of  random 


2The  dimensional  variables,  denoted  with  a  hat,  have  been  non-dimensionalized  using: 

t  :  t  A  /  — 7  :  X  =  xh;  U  =  U  \/ gr  ll :  P  —  pmin  p  (/imax  Pmin) 
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events  u  belonging  to  a  measurable  sample  space  0  and  consider  the  stochas¬ 
tic  velocity,  density  and  pressure  fields:  u(x,  t;uj);  p{x.,t;ui)  and  p(x,  t;oj). 
This  leads  to  stochastic  dynamical  rates-of- change  Cu  and  Cp.  If  these  rate- 
of-changes  are  themselves  uncertain,  for  example  due  to  parameter  or  model 
uncertainties,  then  they  also  depend  explicitly  on  u>.  In  this  study,  we  mostly 
focus  on  uncertainties  arising  due  to  uncertain  initial  conditions.  We  define 
general  stochastic  initial  conditions  as 


u  (x,  0;ca)  =  u0  (x;  co)  ,  xeD,  well, 
p  (x,  0;  co)  =  p0  (x;  co) ,  x  G  V,  co  E  Q, 

and  stochastic  boundary  conditions  as 


U  =  g D  (X,  t]  Co)  , 

<9u  ,  . 

—  =  g^v  (x,t;u;) , 

P  =  9dp  (x,t;a;), 

dP  (  +  \ 

—  gNp  (x,  t\ co) , 


x  G  dVD, 

co  G 

x  G  dVN , 

co  G 

x  G  dVDp, 

co  G 

x  G  dVNp, 

co  G 

(2) 


(3) 


where  the  boundary  conditions  are  separated  into  Dirichlet  and  Neumann 
conditions  for  the  velocity  and  density  fields  (pressure  boundary  conditions 
are  considered  later).  The  resulting  multivariate  stochastic  eqns.  (l)-(3)  de¬ 
fine  the  problem  to  be  solved.  As  in  the  deterministic  case,  specifics  of  the 
solution  depend  on  the  initial  and  boundary  conditions  chosen. 

The  DO  decomposition  of  these  equations  can  be  reconstructed  from 
Sapsis  and  Lermusiaux  [52]  and  appears  in  Sapsis  [51].  For  convenience, 
we  provide  an  abbreviated  summary  of  these  equations  in  Appendix  A. 
In  short,  the  DO  methodology  begins  with  a  generalized  Karhunen-Loeve 
expansion  truncated  to  s(t)  terms,  where  the  number  of  terms  in  the  ex¬ 
pansion  can  generally  be  time-dependent  [53].  The  vector  of  prognostic 
state  variables  3>(x,  t;uj)  =  [u,  p]T  is  decomposed  into  the  sum  of  a  de¬ 
terministic  mean  component  $(x,  t),  with  s  deterministic  modes  $*(x,  t), 
each  mode  multiplied  by  a  stochastic  coefficient  Yi(t]uj).  This  decomposi¬ 
tion  is  first  substituted  into  into  eqns.  (l)-(3).  The  key  DO  condition,  the 
rate-of- change  of  the  stochastic  subspace  is  dynamically  orthogonal  to  the 
subspace  itself,  is  then  utilized.  Orthogonality  is  defined  by  the  spatial  inner- 
product  (a,  b)^  =  fDY2i(atbi)  f°r  arbitrary  vectors  of  spatial  functions 
a  =  [a1,  a2, . .  .]T  and  b  =  [61,  b2, . .  ]7 .  In  general,  we  note  that  this  definition 
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of  the  inner  product  assumes  that  the  different  components  of  the  vector  (in 
our  case,  the  state  vector)  have  been  properly  normalized  [35,  54],  This  is 
not  guaranteed  from  the  simple  deterministic  non-dimensionalization  used  in 
eqn.  (1).  In  fact,  an  additional  stochastic  normalization  is  usually  needed, 
reflecting  the  stochastic  initial  and  boundary  conditions.  After  some  math¬ 
ematical  manipulation  (see  Appendix  A  and  use:  [52,  51])  time-evolution 
equations  for  the  mean,  modes,  and  stochastic  coefficients,  which  are  com¬ 
pletely  determined  by  the  dynamics,  are  obtained.  A  major  contribution  of 
this  manuscript  is  to  derive  efficient  discretizations  in  time  and  space  for  these 
equations  and  to  evaluate  the  resulting  computational  schemes  through  a  set 
of  new  benchmarks  and  applications  for  stochastic  Boussinesq  dynamics. 

3.  Semi-implicit  Time  Discretization 

Solving  the  deterministic  version  of  the  system  of  equations  (1)  implicitly 
in  time  often  requires  not  only  a  large  matrix  inversion  at  each  time-step, 
but  also  iterations  at  each  time-step  to  deal  with  the  non-linear  advection 
terms,  e.g.  [13].  Discretizing  their  stochastic  version  (l)-(3)  using  a  brute- 
force  Monte-Carlo  scheme  would  have  similar  costs  per  realizations,  hence  a 
total  cost  equal  to  that  of  the  deterministic  version  but  multiplied  by  the  size 
of  the  ensemble.  If  a  DO  decomposition  is  used,  solving  the  DO  system  (A.  5)- 
(A.12)  implicitly  would  require  a  matrix  inversion  (s2  +  s  +  1)  times  larger 
than  for  (1)  since  the  mean  and  the  modes  are  coupled  through  the  pressure 
and  non-linear  advection  terms  (the  number  of  pressure  equations  are:  s 2  for 
p/ s,  s  for  pi  s  and  1  for  p).  While  it  is  possible  to  solve  such  systems,  our  goal 
here  is  to  discretize  (A.5)-(A.12)  such  that  the  equations  decouple,  resulting 
in  a  simpler  and  efficient  solution  scheme.  This  section  describes  how  this 
decoupling  is  achieved.  First  we  explain  why  we  treat  some  terms  explicitly 
and  others  implicitly.  We  then  define  new  pseudo-stochastic  pressures  that 
substantially  reduce  computational  costs,  develop  Projection  methods  for 
DO  equations  so  as  to  split  the  velocity  and  pressure  terms,  and  present 
time  marching  schemes  for  the  stochastic  coefficients.  The  complete  time 
discretizations,  combining  these  time  marching  schemes  with  the  projection 
schemes,  are  summarized  at  the  end.  The  spatial  discretizations  of  physical 
space  and  of  the  stochastic  subspace  are  given  in  Sect.  §4. 


3.1.  Explicitly  and  Implicitly  Treated  Terms 

Much  of  the  decoupling  can  be  achieved  by  treating  some  terms  explicitly, 
resulting  in  a  semi-implicit  scheme.  First,  we  choose  to  advance  the  stochas¬ 
tic  coefficients  explicitly,  because  then  Cy^.  and  M y,ymyn  can  be  treated  as 
constants  when  evolving  the  mean  and  the  modes,  and  no  iteration  will  be  re¬ 
quired  between  solving  (A. 5),  (A. 8)  and  (A. 11).  Somewhat  similarly,  we  treat 
the  inner  product  terms  (Q,,  &j)v  &j  in  (A. 8)  explicitly  to  avoid  iterations. 
Next,  we  treat  the  non-linear  advection  terms  explicitly,  which  is  often  done 
in  the  Projection  method  community  (for  e.g.  Guermond  et  al  [19]).  This 
does  impose  a  stability  constraint  on  the  time-step  size,  which  is  now  limited 
by  the  Courant  -  Friedrichs  -  Lewy  (CFL)  condition.  Third,  we  treat  the  lin¬ 
ear  diffusion  terms  implicitly  because  they  do  not  couple  the  equations  and 
the  resulting  diagonal-dominant  matrices  can  be  inverted  efficiently.  While 
they  could  also  be  treated  explicitly,  this  imposes  a  much  harsher  stabil¬ 
ity  constraint  that  could  result  in  very  small  time  steps.  Thus,  to  partially 
decouple  the  evolution  equations,  we  advance  the  stochastic  coefficients,  in¬ 
ner  product  terms  and  non-linear  advection  explicitly.  Unfortunately,  these 
equations  are  still  coupled  through  the  pressure,  which  is  discussed  in  the 
next  section. 

3.2.  The  Pseudo-Stochastic  Pressures 

In  this  section  we  first  review  the  explicit  treatment  of  pressure  which 
results  in  a  scheme  requiring  s2  +  s+l  solutions  of  stochastic  Pressure  Poisson 
equations  (PPEs)  per  time  step.  Then,  we  discuss  our  new  definition  of 
pseudo-stochastic  pressures  that  reduces  the  expense  to  s  +  1  and  show  that 
it  is  a  valid  definition  that  does  not  change  the  order  of  accuracy. 

One  approach  for  handling  the  pressure,  which  was  adopted  in  Sapsis  and 
Lermusiaux  [52] ,  is  to  treat  it  explicitly.  This  approach  takes  advantage  of  the 
fact  that  the  full  stochastic  pressure  can  be  recovered  at  any  time  instant  by 
taking  the  divergence  of  (A. 5)  and  (A. 8),  inserting  the  decomposition  (A. 2), 
and  using  the  divergence-free  form  on  continuity,  noting  that  V  •  ^  =  0.  The 
result  gives  the  stochastic  Pressure  Poisson  equations  (PPEs)  for  our  system 

V2p  =  -V  •  [V  •  (flu)  -  pe9}, 

VVi  =  -V  •  [V  •  (uuj)  +  V  •  (uju)  -  pte9},  (4) 

X?2pi:j  =  -V  •  [V  •  (UjUi)  +  V  ■  (iiiUj)]. 
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A  disadvantage  of  this  explicit  pressure  approach  is  that  it  is  expensive;  to 
recover  the  full  stochastic  pressure,  1  +  s  +  s2  Poisson  equations  need  to  be 
inverted,  and  this  would  often  be  the  dominating  cost  of  the  scheme.  Oceanic 
applications  are  expected  to  require  s  r 0(  102  —  103)  [35],  which  would  be 
very  expensive.  Another  disadvantage  is  that  the  velocity  computed  with  an 
explicit  scheme  will  not  be  divergence-free  after  each  time  step. 

We  can  reduce  the  number  of  PPEs  to  s  +  1  by  defining  new  pseudo¬ 
stochastic  pressures.  The  purpose  of  the  pressure  in  divergence- free  flows 
is  to  enforce  continuity.  In  our  stochastic  equations,  continuity  needs  to  be 
satisfied  by  the  mean  and  each  modal  velocity  field  independently  (we  assume 
that  the  divergence-free  continuity  equation  is  exact,  without  any  errors  in 
its  form).  Also,  each  of  these  velocity  fields  only  needs  a  single  scalar  held  in 
order  to  satisfy  the  continuity  constraint.  By  inspection  of  equations  (A. 5) 
and  (A. 8),  we  therefore  define  new  pseudo-stochastic  pressures,  which  are  a 
combination  of  the  mean,  linear-,  and  quadratic-modal  pressures: 

P 

Pi 

With  this  definition,  the  quadratic  modal  pressures  are  eliminated  from  (A. 5) 
and  (A. 8).  Thus,  to  evolve  the  mean  and  modes,  we  no  longer  need  to  solve 
for  the  quadratic  pressures.  However,  substituting  (5)  into  the  equation  for 
the  evolution  of  the  stochastic  coefficients  (A.  11),  we  fold  that  the  second 
term  on  the  right-hand-side  of  (A.  11), 

(V Prrm  T  V  •  (unUm),  U (  YmYn  Cymyn) 

retains  the  projection  of  the  quadratic  stochastic  pressure  terms  in  the  sub¬ 
space.  At  first,  this  would  indicate  that  the  quadratic  modal  pressures  are 
still  needed,  but  for  commonly  used  boundary  conditions,  the  projection  can¬ 
cels,  i.e.  the  inner  product  (Vprm,  Ui)v  is  zero.  We  refer  to  Sapsis  et  al  [54] 
for  the  proof.  Because  of  this  property,  the  quadratic  stochastic  pressure 
term  in  (A.  11)  can  be  dropped  without  any  penalty.  Thus,  by  defining  new 
pseudo-stochastic  pressures  (5),  we  have  shown  that  we  reduced  the  number 
of  PPEs  from  s2  +  s  +  1  to  the  expected  s  +  1. 

3.3.  Projection  Methods  for  the  Mean  and  Modes 

To  obtain  a  numerically  divergence-free  velocity,  we  use  a  Projection 
method.  A  large  number  of  different  Projection  methods  exist;  for  a  re¬ 
cent  review,  see  [19].  Projection  methods  are  known  for  excellent  efficiency, 


=  P  +  Cy,y:i  Pij, 

=  pi  +  Cy^y.  M YjYmYn  Prm- 
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but  the  proper  specification  of  boundary  conditions  remains  a  long-standing 
issue.  While  advances  in  this  area  are  still  being  made  [16,  55],  we  have  cho¬ 
sen  to  use  the  “incremental  pressure-correction  scheme  in  rotational  form” 
proposed  by  Timmermans  et  al  [60],  which  has  a  proven  temporal  accuracy 
[19].  We  first  summarize  the  classic  versions  of  the  scheme,  then  adapt  them 
for  the  mean  and  modes. 

Classic  Projection  Scheme.  Intermediate  velocities  are  first  solved  for 
using  a  first  or  second  order  time-accurate  discretization.  In  both  cases,  the 
contributions  from  known  velocities  (at  previous  or  intermediate  time  steps) 
are  written  as  F(ufc*,pfc*),  where  the  time  instant  tk*  determines  the  order 
of  the  scheme  in  time.  Specifically,  for  first  order  in  time,  k*  =  k  —  1,  and 


^  -  ^=W  =  _Vpt"1  +  x  e  v> 

ufc  =  gD,  x  G  &Dd , 
dnk 

— —  =  gjv,  x  G  dVN. 
on 

This  is  followed  by  the  computation  of  the  pressure  correction,  0,  so  as  to 
satisfy  continuity, 


V2d  = 


V  •  uk 

At 


dd 


=  0,  x  G  &Dd , 


dn 

9  —  9pi  x  £  &DN , 


where  gp  is  a  the  prescribed  Pressure  difference  at  the  boundary.  With  this 
correction,  the  full  pressure  and  velocity  fields  can  be  recovered  at  the  next 
time  step,  i.e. 

uk  =  uk  -  V9k, 

pk  =  pk-]  +  ek  -  uv  ■  uk. 


The  first-order  time-accurate  discretization,  k *  =  k  —  1,  uses  the  con¬ 
tributions  from  known  velocities  (at  old  time  steps)  as  F(ufc-1,  pk~l).  The 
second-order  time  integration  scheme  is  obtained  by  appropriately  choosing 
the  guess  values  (»)fc*  with  tk-i  <  tk*  <  tk,  and  using  a  higher  order  back¬ 
wards  differencing  method.  For  additional  properties  of  such  schemes,  we 
refer  to  Timmermans  et  al  [60],  Guermond  et  al  [19]. 
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Projection  Scheme  for  the  Mean.  We  evolve  the  mean  fields  modifying 
the  classic  projection  method  to  account  for  the  moments  of  the  stochastic 
coefficients  and  of  the  chosen  explicit  and  implicit  terms  (Sect.  §3.1).  Starting 
from  the  PDEs  (A. 5)  for  the  mean,  we  obtain: 


afc 


2  =  fc 


- - ^=V  u 

A  t  jGf 


Vz9 


2  nk 


pk 


A  t  ScV(fr 


ufc 

f 


V2pk 


^  -  {v  •  (uu)}fc*  -  vr +pfev 

-  C$iYj  {V  •  (u,u4)}fc* , 

-J-V  • 

At 

-  A tW9k, 

f- 1  +  9k  -  nW  •  ufc, 

-  {V  •  (w)r  -  {V  •  (u^)}fc* , 


(6a) 

(6b) 

(6c) 

(6d) 

(6e) 


with  deterministic  boundary  conditions: 


89 

x  G  8VD, 

On? 

to 

?|: 

o 

d 

x  e  dVN , 

dn 

=  g  n,9  =  gP , 

Pk  =  g  Dp, 

x  G  dVDp 

8pk 

r\  &Np  5 

an 

x  G  dVNp 

As  above,  k*  =  k  —  1  gives  the  first  order  time-accurate  scheme,  and  for 
second  order  accuracy  we  refer  to  Guermond  et  al  [19].  We  note  that  a  main 
difference  between  the  classic  scheme  and  the  above  new  DO  mean  scheme 
are  the  presence  of  the  covariances  and  third  moments  of  the  coefficients 
Yf  s.  When  to  evaluate  these  is  discussed  in  Sect.  §3.4  and  the  order  in  which 
the  coupled  differential  equations  for  the  mean,  mode  and  coefficients  are 
integrated  is  presented  in  Sect.  §3.5. 

Projection  Scheme  for  the  Modes.  As  for  the  mean,  the  modes  are  evolved 
by  modifying  the  classic  projection  method  for  the  modal  PDEs  (A. 8).  We 
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obtain: 


u? 


A  t  ^ 


V 2u? 


u 


k- 1 


At 


{V  ■  (UiQ)}*'  -  {V  ■  (flUj)}*'  -  Vjf  +  pfe« 


^YiYj  M-YjYmYn  '  (UnW-m)} 


i-l,  k* 

'YiYj 

(Qn^uf, 


r  k* 

lY 

,fe*  _  ,fc* 


V^i  =^V  •  uf , 

U?  =  U?  -  AfWf, 

$=Afe“1+0f-i'V-u£, 


p\ 


At  5c^/G^ 


vVf  =  p: 


fc-i 


At 


{V  ■  (uj/i)}fc*  —  {V  •  (up*)}* 


C™  Mtymyn  {V  ■  (u nPm)} 


(Q 


with  boundary  conditions: 


00. 


<9u,-; 

dn 


W  =  g i,D(=  0),  —  =  0,  X  e  dVD , 

=  gi,jv(=  o),  6>i  =  ft)P(=  0),  X  e  dVN, 

pik  =  gi,Dp(=  0),  *edVDp, 

Qn-k 

—  -  »)■  xeMV„, 


(8a) 

(8b) 

(8c) 

(8d) 

(8e) 


(9) 


where,  again,  the  scheme  is  hrst  order  for  k*  =  —  1.  Since  we  focus  on 

the  numerics  of  the  DO  differential  equations,  we  assume  that  the  stochastic 
boundary  forcings  are  null,  i.e.  g^s  are  null  in  eqns.  (9).  We  note  that  prob¬ 
lems  with  stochastic  boundary  forcing  can  always  be  transformed  into  a  prob¬ 
lem  with  deterministic  boundary  conditions  and  suitable  interior  stochastic 
forcing  fields  (see  Chap.  5  of  Sapsis  [51]).  In  the  present  study,  we  do  not 
exemplify  problems  with  stochastic  forcing,  but  we  touch  on  the  numerical 
modifications  required  in  such  cases  (see  §3.4  and  [34]  for  an  example  of 
stochastic  ocean  models).  The  time  evolution  of  the  l^’s  is  discussed  next. 
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3-4-  Time  Integration  Scheme  for  the  Stochastic  Coefficients 

To  integrate  the  Ordinary  Differential  Equations  (ODEs)  (A.  11)  for  the 
stochastic  coefficients,  we  assume  that  all  variables  are  available  at  time  tk-\ 
and  that  we  integrate  forward  to  time  tk- 

First,  we  consider  the  case  where  the  original  governing  differential  equa¬ 
tions  only  contain  uncertain  initial  conditions,  and  no  stochastic  forcing. 
This  case  corresponds  to  the  examples  and  benchmarks  we  consider  later. 
It  allows  the  use  of  classic  time-marching  integration  algorithms  or,  in  other 
words,  appropriate  ODE  solvers.  We  consider  a  given  realization  of  a  coeffi¬ 
cient  Yi  at  time  tk- 1-  To  integrate  to  time  tk,  we  first  define  an  approximation 
to  the  time-rate  of  change  of  this  coefficient  at  a  given  instant  tk- 1  <  t  <  tk 
as  follows: 


dt 


Y(t) 


V  •  (umu)  -  V  •  (uum)  -  Vpm  +  Pm.e9,  u,: 


k~ 

V 


Ymft ) 


V  •  (U mp)  -  V  •  (U pm),  Pi 


Ym(t) 

V 


-  (V  •  (unum),  nf)kv  ( Ym(t )  Yn(t )  -  C*Tyn) 

-  (V  •  (u npum),pf)%  (Vm(t)  Yn(t)  -  C^myn)  , 


(10) 


where  (. )k  indicates  that  the  modal  quantities  are  estimates  of  their  values 
at  time  t.  The  numerically  exact  option  is  to  choose  tk-  =  t,  while  the 
cheapest  is  to  take  tk-  =  tk- 1  since  at  initial  time,  all  modal  quantities  are 
available  from  the  previous  time-step.  Using  this  derivative  definition,  we 
have  employed  several  time-marching  schemes  of  varying  order  to  advance 
the  Yi  in  time,  including:  a  low-storage  4th-order-accurate  explicit  Runge- 
Kutta  integrator  of  the  form 


y(°)  =  Y(t),  y{0)  =  0, 

y(m)  =  amy^  + At 

dt 


y(m-!)  ,t+cmAt 


y(m)  _  y(m- 1)  +  ^(m) 

Y(t  +  At)  =  y(5), 


for  me  [1 ...  5] , 
for  me  [1 ...  5] , 


where  the  coefficients  am,bm,cm  are  given  in  Carpenter  and  Kennedy  [3];  a 
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2nd-accurate  explicit  Runge-Kutta  scheme  (Heun’s  version) 


y(0)  =  Y(t)  +  At 


dY 

dt 


Y(t  +  At)  =  Y{t)  + 


'  Y(t) 

At  I  dY  I 


dt 


+ 


Y(t) 


dY 

dt 


y(  o) 


and  the  first-order-accurate  explicit  Euler  scheme 


dY 

Y(t  +  At)  =Y(t)+At  — 


Yt 


For  the  Euler  scheme,  k~  =  k  —  1.  For  each  stage  of  the  two  RK  schemes, 
'-A-  is  evaluated  using  (10).  If  tk-  =  t,  the  modal  quantities  (the  rates  for  the 
Yds)  are  advanced  to  intermediate  times  using  the  modes  and  mean  PDEs 
above  (which  is  expensive).  If  k~  =  k  —  1,  the  mean  and  modal  quantities 
are  not  advanced,  but  the  Ym(t)'s  are  still  updated  at  intermediate  times.  Of 
course,  the  formal  order  of  accuracy  of  that  scheme  is  limited  by  these  rates 
kept  constant  as  will  be  shown  later  (§6).  While  C YmYn  could  be  recalculated 
at  each  time  level,  our  simulations  showed  better  results  if  they  were  kept  at 
the  same  time  k~  as  modal  quantities  (in  that  case,  coefficients  and  subspace 
remain  consistent). 

Second,  we  consider  the  case  where  the  original  governing  equations  con¬ 
tain  stochastic  forcing  in  the  form  of  zero-mean  Wiener  processes,  added 
linearly  to  the  deterministic  form  of  the  PDEs  (1).  These  stochastic  forcings 
then  appear  in  the  governing  equation  (A. 11)  for  the  stochastic  coefficients. 
Without  entering  details,  the  above  marching  schemes  are  then  augmented 
with  stochastic  integrals  over  tk  —  tk_ 1  for  the  contribution  of  these  forcings 
[26,  22],  If  the  stochastic  forcings  are  of  multiplicative  form  in  the  original 
governing  equations,  then  expectations  and  coupled  integrals  of  stochastic 
terms  are  in  general  also  in  the  mean  and  mode  equations  of  Sect.  §3.3.  This 
is  not  considered  in  the  present  study. 


3.5.  Complete  Time  Integration  Scheme 

We  now  summarize  the  complete  time-discretization  scheme  from  tk-i  to 
tk.  Since  we  have  decoupled  the  equations  (6),  (8),  and  (10),  the  order  in 
which  they  are  solved  is  not  important.  In  fact,  (6),  (8),  and  (10)  could  be 
solved  in  parallel.  Presently,  we  employ  the  following,  serial  approach: 
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1.  Calculate/extrapolate  the  statistics  (Cymy„,  M^^yJ  to  the  approxi¬ 
mated  times  k*,  k~  and  store  for  later  use 

2.  Calculate/extrapolate  the  advection  terms  (Vuu,  VuUj+VujU,  Vu;Uj) 
to  the  approximated  times  k*,  k~  and  store  for  later  use 

3.  Advance  the  Yd s  using  (10),  and  one  of  the  ODE  solvers  in  §3.4 

4.  Advance  the  mean  u  using  (6) 

5.  Advance  the  modes  tp  using  (8). 

For  the  modes  and  mean,  we  choose  k*  =  k~  =  k  —  1,  resulting  in  a  first 
order  accurate  scheme  for  the  present  applications  (§6).  For  second  order 
accuracy,  we  would  use  k *  as  defined  in  Timmermans  et  al  [60],  Guermond 
et  al  [19].  For  the  stochastic  coefficients,  a  higher  order  ODE  solver  may  be 
used  in  step  3  to  reduce  the  magnitude  of  the  integration  error,  but  the  order 
of  accuracy  of  that  step  is  still  limited  by  the  choice  of  k~,  and  hence  will  be 
first  order  accurate. 

4.  Spatial  Physical  and  Stochastic  Subspace  Discretizations 

In  this  section,  we  start  by  describing  the  2D  spatial  discretization  of 
(6)  and  (8)  on  a  structured  grid  (§4.1).  We  employ  a  standard  conservative 
finite  volume  discretization  of  second-order.  Special  treatment  is  needed  for 
the  advection  by  the  modes  since  they  are  basis  vectors  in  the  stochastic 
subspace  and  thus  do  not  have  a  preferential  direction.  We  finally  discuss 
the  discretization  of  the  s-dimensional  probability  subspace  in  §4.2. 

4-1.  Spatial  Discretization  of  the  Physical  Space 

The  domain  V  is  discretized  into  a  finite  number  of  non-overlapping  con¬ 
trol  volumes.  Presently,  we  use  rectangular  control  volumes,  which  form  a 
structured  Cartesian  grid  that  has  uniform  spacing  in  the  x  and  y  directions. 
Different  choices  exist  for  the  relative  placement  of  velocity  and  pressure  con¬ 
trol  volumes.  Here  we  choose  to  use  a  standard  staggered  C-grid  [13,  42], 
where  the  u-  and  u-velocity  control  volumes  are  displaced  half  a  grid-cell  in 
the  x-  and  ^/-directions  relative  to  the  pressure  and  density  control  volumes, 
respectively. 

4-1.1.  Diffusion  Operator 

The  diffusion  operator  V2  is  discretized  by  using  central  boundary  fluxes. 
That  is,  for  the  held  <f  at  the  center  (x,  y )  of  a  control  volume  boundary,  we 
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have  the  following  boundary  flux  in  the  ^-direction  (similar  for  ^/-direction) 

50  =  0(s  +  j£jy) 

9X  (*,y)  Aa: 

which,  on  a  structured  grid,  is  exactly  equivalent  to  the  second-order  accurate 
central  finite-difference  scheme.  For  the  advection  operator,  a  simple  second- 
order  central  flux  is  however  well-known  to  be  unstable  (e.g.  [4]),  and  needs 
more  careful  treatment,  as  described  next. 


4-1.2.  Advection  Operator 

For  advection  by  a  velocity  component  u,  we  use  a  standard  Total  Vari¬ 
ation  Diminishing  (TVD)  scheme,  with  an  monotonized  central  (MC)  sym¬ 
metric  flux  limiter  [61].  The  scheme  can  be  written  for  a  variable  r)  as: 


F(Vi-±)  =  ut_i 


Vi  +  rjj-i 
2 


(11) 


where  the  MC  slope  limiter  T(r)  is  defined  as 


T(r)  =  max  <  0,  min  min 


1  +  r 


and  the  variable  r  as 


i  = 
1  2 


l  (ut_  1  +  IV 1  )  (Vi- 1  +  Vi- 2)  +  I  (ui-l  -  )  (Vi+1  +  Vi) 

u.i_i(vi  -Vi-0 


where  ut_i  is  used  without  interpolation  for  the  density  advection  while  a 
second-orcfer  central  scheme  is  used  for  the  non-linear  u  and  v  advection.  For 
more  information  about  TVD  schemes,  we  refer  to  [39]. 

A  possible  issue  with  using  this  scheme  for  DO  equations  arises  from  the 
realization  that  the  absolute  value  of  the  velocity  |u|,  is  a  function  of  the 
full  velocity  u  +  Ylul  which,  depending  on  the  specific  realization,  may  be 
either  positive  or  negative;  in  other  words,  |u|  is  positive,  but  stochastic. 
Fortunately,  we  never  need  the  full  velocity  to  evolve  the  mean  and  modes 
in  (A. 5)  and  (A. 8)  (see  also  (6)  and  (8)).  In  fact,  in  the  case  of  the  mean 
velocity  u,  its  absolute  value  is  deterministic.  Therefore,  the  advection  of 
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the  mean  velocity  by  the  mean  velocity  u  •  Vu,  and  the  advection  of  the 
velocity  modes  by  the  mean  velocity  u- Vu*  can  use  the  classic  TVD  method 
without  modification.  Advection  of  the  mean  by  the  modes  u,  •  Vu  and  of 
the  modes  by  the  modes  Uj  •  Vu.,-,  however,  need  additional  consideration. 
Similar  statements  apply  for  the  advection  of  the  mean  density  and  of  the 
density  modes  by  either  the  mean  velocity  (classic  scheme  is  fine)  or  by  the 
modes  (additional  considerations  are  needed). 

Here  we  propose  three  arguments  for  three  different  advection  schemes 
that  can  be  used  for  these  “advection  by  the  modes”  terms.  First,  if  we 
examine  the  equations  from  the  perspective  of  the  numerical  scheme  only,  a 
preferential  advection  direction  will  be  present.  In  this  case  we  simply  use 
the  TVD  scheme  unmodified.  Next,  we  argue  that,  since  the  stochastic  coeffi¬ 
cients  are  zero  mean,  then  the  probability  of  Yj  <  0  is  equal  to  the  probability 
that  Yj  >  0.  This  suggests  that  u*  should  not  have  a  preferential  direction 
of  propagation,  in  which  case  a  central  differencing  advection  scheme  (CDS) 
could  be  used.  Last,  recognizing  that  the  CDS  scheme  may  cause  oscillations, 
we  still  wish  to  limit  the  flux  in  some  way.  A  direct  approach,  then,  is  to 
use  the  TVD  scheme  in  both  directions,  and  average  the  results.  That  is, 
the  present  sign  of  the  modal  velocity  is  used  first  to  calculate  the  advective 
terms,  then  the  negative  of  the  modal  velocity  is  used,  and  the  two  results 
are  averaged.  We  call  this  the  symmetric  TVD  or  TVD*  scheme.  Thus,  we 
have  three  potential  schemes  for  advection  by  the  modal  velocities. 

The  three  proposed  schemes  are  tested  in  §6.2,  where  we  find  that  the 
TVD*  scheme  performs  best.  Improvement  on  this  scheme  is  still  possible, 
since  we  also  find  that  very  small  oscillations  can  remain.  A  proper  flux- 
limited  advection  scheme  may  be  derived  by  looking  at  the  characteristics 
of  the  hyperbolic  parts  of  the  full  system.  However,  since  the  system  has 
(s  T  1)  x  d  equations,  where  d  is  the  dimension  of  the  problem,  this  analysis 
is  left  for  future  research. 

4-2.  Discretization  of  Stochastic  Subspace 

The  stochastic  coefficients  exist  in  an  .s- dimensional  space,  which  can  be 
potentially  large,  from  0(10)  to  O(103)  based  on  our  experience.  In  most 
cases,  there  is  also  no  strict  bound  on  the  value  that  a  stochastic  coeffi¬ 
cient  can  take.  Thus,  evenly  dividing  the  ,s- dimensional  space  is  not  feasible. 
To  discretize  the  uncertainty  subspace,  other  schemes  are  thus  used.  They 
include:  (a)  non-uniform  discretizations  of  the  subspace,  either  using  struc¬ 
tured  or  unstructured  grids,  possibly  using  schemes  based  on  finite-volumes 
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or  finite-elements  [18,  43];  (b)  solve  a  discretized  version  of  the  PDEs  for  the 
probability  densities  of  the  coupled  s  coefficients,  e.g.  solve  Fokker-Planck 
equations  [52];  (c)  parameterize  the  probability  space,  either  using  Polyno¬ 
mial  Chaos  [18,  69,  44,  15],  in  our  case  extended  to  time-dependent  polyno¬ 
mials  [e.g.  17],  or  other  parameterizations  such  Gaussian  mixtures  [46,  58]  or 
particle  filters  [9],  and  (d),  use  a  Monte-Carlo  approach  [37,  32], 

We  have  employed  a  few  of  these  schemes.  In  the  present  manuscript, 
we  only  illustrate  the  use  of  a  Monte-Carlo  scheme  using  for  each  realization 
the  time-integration  described  in  Sect.  §3.4.  In  general,  the  expected  error 
for  the  mean  and  covariance  is  of  O  ^-4=^,  where  q  is  the  number  of  samples. 
For  efficient  results,  it  is  important  that  samples  are  generated  in  regions 
where  the  probability  is  relatively  high  and  importance  sampling  [9]  can  be 
employed.  At  initial  time  to,  the  distribution  of  the  Yj  is  generated  using  the 
specified  initial  probability  density  given  by  eqns.  (2).  Here,  our  focus  is  on 
numerical  schemes  and  we  restrict  ourselves  to  simple  distributions  such  as 
Gaussians  or  Dirac  functions  for  numerical  tests. 

Thus,  we  discretize  Y[  by  generating  q  samples,  and  forming  a  qxs  dimen¬ 
sional  matrix  Yri.  Each  row  in  Yr ^  corresponds  to  one  of  the  samples,  and 
each  column  corresponds  to  one  of  the  modes.  During  the  time-integration 
step,  each  sample  of  Yr^  is  then  advanced  using  (10),  which  is  done  efficiently. 
In  all  cases  computed  so  far  q,  can  be  large,  e.g.~  (9(104  — 105),  but  still  suffi¬ 
ciently  small  such  that  advancing  (10)  for  every  sample  is  not  the  dominating 
cost  of  the  whole  scheme. 

A  drawback  to  this  Monte-Carlo  approach  is  that  rare  events  will  not 
be  captured  unless  a  very  large  number  of  samples  are  used.  Alternative 
methods,  such  as  Gaussian  mixtures  [57,  58]  or  other  approaches  mentioned 
above  can  then  be  used  as  an  alternative. 

5.  Implementation  Details 

In  this  section  we  describe  selected  implementation  details  that  arise  due 
to  numerical  issues.  In  particular,  we  discuss  how  to  deal  with  possibly 
poorly  conditioned  covariance  matrices  in  the  stochastic  subspace,  as  well 
as  the  orthonormalization  of  the  modes  and  decorrelation  of  the  stochastic 
coefficients. 
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5.1.  Dealing  with  a  Singular  Covariance  Matrix 

The  covariance  matrix  may  be  singular  or  poorly  conditioned  if  one  or 
more  of  the  stochastic  coefficients  have  zero  or  very  small  variance  com¬ 
pared  to  other  modes.  This  situation  can  arise  naturally  if  a  problem  has 
deterministic  initial  conditions,  with  uncertainty  introduced  through  forcing, 
boundaries,  parameters,  numerical  uncertainties,  or  other  causes.  In  this  ex¬ 
ample,  the  initial  covariance  matrix  is  simply  zero,  for  which  the  inverse  is 
not  defined.  Special  treatment  is  thus  needed  for  such  cases  since  the  inverse 
of  the  covariance  matrix  is  required  in  (8). 

Fortunately,  this  problem  is  also  common  in  data  assimilation  [1,  37] 
and  is  relatively  easily  resolved  using  generalized  Moore-Penrose  inversions. 
First,  in  the  particular  case  of  the  DO  equations,  the  inverse  of  the  covariance 
matrix  is  multiplied  by  the  third  moments  in  (8), 


l,k 


Mf 


(12) 


which,  for  most  physical  processes,  goes  to  zero  for  the  eigenvalues  of  Cy.Y- 
that  go  to  zero.  To  obtain  a  numerically  stable  estimate  of  the  inverse  of  the 
covariance  matrix,  we  simply  perform  a  generalized  or  pseudo  inverse,  setting 
the  singular  values  less  than  a  defined  tolerance  equal  to  zero  or  to  that  given 
tolerance.  This  results  in  a  stable  numerical  simulation,  as  exemplified  for 
the  lock-exchange  validation  benchmark  in  §6.1. 


5.2.  Orthonormalization 

The  DO  equations  require  orthonormal  modes,  therefore  it  is  important  to 
maintain  this  property  numerically  when  integrating  the  equations  over  time. 
Starting  with  initially  orthonormal  modes,  orthonormality  is  maintained  over 
time  when  applying  the  DO  condition  in  continuous  mathematics,  since 


3/V 


dt 


d<i\ 

dt 


$ 


v 


<I> 


?*3 

dt 


=  0, 


v 


where  the  DO  condition  =  Vi,j  G  [l,2,...,s]  was  used. 

At  the  discrete  level,  this  property  is  maintained,  but  only  in  a  weak  sense, 
up  to  truncation  and  round-off  errors.  Therefore,  even  if  the  modes  are 
orthonormal  at  a  given  time-step,  not  all  integration  schemes  over  the  next 
time-step  will  conserve  the  discrete  orthonormality.  To  see  this,  consider  the 
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inner  product  at  time  tk,  integrated  from  time  tk- 1,  at  first  exactly  (without 
any  discretization): 


<*4 ,  =  <**,**>,, , 


=  (  $-'_1  + 


d®, 

dt 


dt ,  fcj"1  + 


d*. 

dt 


dt 


v 


{  1  dt  '  3  dt  I 


=  (S*-1 ,  *J“%  +  A tU 


k-i  ^ 


’  dt 


v 


At(^~\^l)  +  At2'"~3 


dt 


v 


d&i  d<t>, 


dt  ’  dt 


v 


—  Oij , 


where  we  have  used  A t^3-  =  f  Ljjj-d/.  Now,  a  pth-order- accurate  discrete 

approximation  of  will  have  an  associated  error,  ef,  of  0(Atp).  Assuming 
that  the  spatial  inner  product  is  computed  exactly,  this  error  in  the  time 
integration  leads  to  an  error,  £|,  in  (<frj ,  $j)^dlscrete  =  (<3?* ,  (E»J).p  +  £|  at 
time  k.  To  find  the  magnitude  of  this  time  discretization,  £k,  we  assume 
to  start  with  an  exactly  orthonormal  inner  product  at  tk-i,  then  by  discrete 
integration  to  tk  we  have: 


0  numerically  using  DO  condition 


(#, ,  +£*  =  (#J_1 ,  «*-%  +Ai  (  a>‘-‘ ,  ^  +  ej 


P 


0  numerically  using  DO  condition 


At  (  ^ 3  +  ek 

1  ’  <9t  2 


A  2  /  d*i 

At<~ r 


k  d®L 

+  £*,—  +  £ 


dt 


v 
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0  because  exact 


=  **i\  +At2  ( ek  **i\ 

•  Af\dt'dt/V+At\£^dt/V 


At2  /£fc 
\  1  ’  dt 


+  <4  >  > 


2Af2e>(Atp)e>(Af~1)  +  At20(At2p), 
0(Atp+1). 


Therefore,  unless  the  discrete  time  rates  of  change  are  also  made  orthogonal, 
or  the  sum  of  the  three  terms  above  are  discretely  equal  to  zero  (assuming  the 
first  one  is  zero),  the  error  in  the  orthonormality  will  always  be  At  smaller 
than  that  of  the  numerical  scheme. 

The  orthonormality  can  be  corrected  either  directly  by  enforcing  that 
)  yfy  +  £i^  =0  numerically  during  time  integration,  or  indirectly 

by  enforcing  that  (<£>,; ,  =  8$  numerically  at  the  end  of  a  time 

step.  In  either  case,  it  is  important  to  remember  that  the  summation  Yr^<&, 
produces  specific  realizations,  and  changing  the  basis  without  modifying  the 
coefficients  will  change  the  specific  realizations.  Because  T,  and  Yri  are 
linked,  various  schemes  for  performing  the  re-orthonormalization  exist.  These 
schemes  are  briefly  discussed  in  Appendix  B. 


6.  Numerical  Applications 

In  this  section  we  present  four  benchmarks  used  to  verify  and  study  the 
implementation  described  above.  To  ensure  that  the  numerical  implemen¬ 
tation  is  solving  the  desired  equations,  we  compare  the  stochastic  code  to  a 
deterministic  code  (§6.1)  for  a  version  of  the  lock-exchange  dynamical  prob¬ 
lem  [21],  Next,  we  examine  the  different  advection  schemes  proposed  for  DO 
equations  in  §4.1.2,  using  a  symmetric  version  of  the  lock-exchange  problem 
(§6.2).  Then,  we  study  the  spatial  and  temporal  convergence  using  the  lid- 
driven  cavity  flow  (§6.3).  Finally,  we  study  the  discretization  of  the  stochastic 
coefficients  using  the  flow  over  a  square  cylinder  in  a  confined  channel  (§6.4). 

For  evaluating  the  errors  in  the  benchmarks,  we  use  the  L2  norm.  At 
a  single  time  instance,  the  norm  is  ||$||2  =  \J (3b  $ )v  for  the  whole  state 

vector,  or  ||0||2  =  \f  Jv  (j)2dD  for  a  single  component.  For  the  convergence 
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x=  (-1,1)-, 


un  =  0,  du/dn  =  0 


Light  Heavy 


x=  (-1,-1)  x  =  (1,-1) 

Figure  1:  An  initial  barrier  separating  light  and  heavy  fluid  is  removed,  and  the  flow 
is  allowed  to  evolve.  Uncertainty  in  our  studies  originates  from  not  knowing  the  initial 
density  differences  between  the  fluids.  This  benchmark  is  used  to  verify  the  correctness  of 
the  implementation,  not  the  DO  methodology. 


studies  (§6.3)  we  also  integrate  over  time,  using  ||3>||^  =  y  J QT  ($,  <&)v  dt  for 
the  state  vector,  and  ||Y)||2  =  \J E  [KWi]  dt  for  the  stochastic  coefficients. 

6.1.  Lock-Exchange  Verification  Benchmark 

The  purpose  of  this  benchmark  is  to  verify  the  numerical  implementation. 
For  example,  in  theory,  the  definition  of  the  stochastic  pseudo-pressure  does 
not  introduce  any  additional  error  in  the  projection  method.  However,  it  is 
necessary  to  verify  that  its  numerical  implementation  is  sufficiently  accurate, 
and  that  it  is  still  solving  the  correct  equations.  Ideally,  problems  with 
analytical  solutions  should  be  used  to  verify  a  code,  however,  constructing  a 
valid  analytical  solution  for  evolving  (6), (8),  and  (10)  with  multiple  stochastic 
modes  and  coefficients  is  neither  trivial,  nor  does  it  lend  itself  to  compact 
expressions.  In  this  section  we  address  this  problem  by  defining  a  numerical 
benchmark  that  can  be  used  to  verify  the  numerical  implementation.  We 
then  run  this  benchmark  for  our  implementation. 

6.1.1.  Lock- Exchange  Verification  Benchmark:  Setup 

We  verify  our  stochastic  DO  code  by  comparing  it  to  a  deterministic  NS 
code  which  has  been  thoroughly  verified  [36].  This  deterministic  code  uses 
the  same  second-order  Finite  volume  scheme  and  first  order  backwards  dif¬ 
ference  Projection  method  as  the  stochastic  code.  While  the  stochastic  code 
is  inherently  more  complicated  with  more  equations,  the  major  differences 
are  in  the  advection  scheme  used  for  advection  by  the  stochastic  modes  (see 
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Figure  2:  Initialization  of  the  lock-exchange  problem  (density  field:  mean,  mode  1  and 
marginal  pdf,  modes  2  and  3).  The  initial  velocity  is  zero,  and  stochasticity  is  introduced 
through  the  first  (of  three)  orthonormal  modes  for  the  density.  The  vertical  length  scale 
used  for  non-dimensionalization  is  the  half-height  of  the  channel  ( h  =  1).  Initializing 
with  four  stochastic  samples  for  the  first  mode,  Y^i  =  [—0.18,  —  0.06,  0.04,  0.20]T,  we 
have  four  possible  initial  conditions  corresponding  to  A p  =  [1,  0.84,  0.74,  0.62].  The  two 
remaining  modes  are  initialized  as  described  in  the  text,  and  do  not  introduce  additional 
stochasticity. 

§4.1.2  and  §6.2),  and  that  we  also  have  to  accurately  evolve  the  stochastic 
coefficients. 

Since  the  deterministic  code  can  be  used  non- intrusively  with  the  Monte- 
Carlo  method  to  generate  an  ensemble  of  realizations,  we  can  compare  re¬ 
alizations  from  the  stochastic  code  directly  to  the  deterministic  realizations. 
To  avoid  generating  a  large  ensemble,  we  prescribe  a  discrete  pdf  for  this 
stochastic  benchmark.  This  benchmark  (Fig.  [1])  is  based  on  the  lock- 
exchange  problem  [21],  where  uncertainty  is  introduced  by  prescribing  four 
possible  density  differences.  That  is,  consider  a  lock  separating  two  miscible 
fluids  of  different  densities.  While  the  difference  between  the  densities  are 
unknown,  it  is  known  that  four  possibilities  of  equal  probability  exists.  An¬ 
other  way  to  look  at  this  benchmark  is  that  we  are,  essentially,  using  a  single 
run  of  the  stochastic  code  to  try  and  replicate  four  different  deterministic 
runs.  While  not  a  practical  use  of  the  DO  method,  it  is  a  challenging  and 
useful  benchmark  problem  for  verifying  its  numerical  implementation. 

To  capture  all  of  the  uncertainty,  three  DO  modes  (s  =  3)  are  needed  for 
the  stochastic  simulation.  Four  density  difference  were  chosen  because  three 
DO  modes  is  the  minimum  number  required  to  have  full  energy  interactions 
between  the  mean  and  modes  of  the  DO  simulation  [52,  51].  Also,  to  fully 
verify  the  implementation,  we  need  to  ensure  that  the  initial  pdf  is  non- 
symmetric  so  that  the  third  moments  are  non- zero  (in  (10)  for  example). 

General  setup:  The  Schmidt  number  is  kept  constant,  Sc  =  1,  and  we 
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present  results  for  two  Grashof  numbers,  Gr  =  1.25  x  106  and  Gr  =  4  x  104, 
(although  more  were  studied).  The  four  possible  density  differences  of  equal 
probability  are  A p  =  [1,  0.84,  0.74,  0.62],  with  the  initial  density  profile 
prescribed  by 

P(x,y,t  =  0)  =  -^tanh(2  x/lp), 

where  we  take  lp  =  1/64.  Initially  the  velocity  is  zero  everywhere,  and  free- 
slip  boundary  conditions  are  used  at  the  boundaries  of  the  domain  (Fig.  [1]). 
The  domain  is  discretized  using  Ax  =  Ay  =  1/256,  and  At  =  1/512,  which 
is  sufficient  resolution  for  these  Grashof  numbers  [21]. 

Mean  initialization:  The  mean  density  profile  (Fig.  [2]-a)  uses  the  hyper¬ 
bolic  tan  function  specified  above  with  the  mean  density  difference  A p  =  0.8. 
The  mean  pressure  and  velocity  are  zero  everywhere,  initially. 

Mode  initialization:  The  density  profile  for  the  first  mode  is  the  normal¬ 
ized  hyperbolic  tan  profile  above  (Fig.  ]2]-b).  The  two  remaining  modes  are 
arbitrary,  since  they  do  not  introduce  initial  uncertainty  (see  below).  They 
are  set  to: 


p2(x,f  =  0) 


p3(x,f  =  0) 


(A p-  |p|)sign(p)|sin(7n/)| 
0 


if  (A p  —  |p|)sign(p)  sin (717/)  >  0 
otherwise, 


(A p-  |p|)sign(p)|  sin(7T7/)| 
0 


if  (A p  —  |p|)sign(p)  sin  (717/)  <  0 
otherwise, 


where  these  are  orthonormalized  numerically  as  described  in  §5.2.  The  pres¬ 
sure  and  velocity  for  all  modes  are  zero  everywhere,  initially. 

Stochastic  coefficient  initialization-.  The  discrete  probability  density  func¬ 
tion  for  the  first  stochastic  coefficient  is  specified  by  using  four  samples 
yr,i  =  [—0.18,  —  0.06,  0.04,  0.20]t,  (Fig.  [2]-c).  The  next  two  modes 

do  not  introduce  additional  uncertainty  because  we  specify  correlated  sam¬ 
ples,  Yrt 2  =  Yr> 3  =  e-  [—0.18,  —0.06,  0.04,  0.20]T,  where  e  is  a  small  constant 
chosen  such  that  JA  Yax(Yri)  =  Var(Kr  l)  numerically.  This  means  that  the 
inverse  of  the  Covariance  matrix  will  be  very  ill-conditioned  (or  numerically 
singular),  so  the  pseudo-inverse  is  required  during  the  initial  stages  of  the 
simulation  (see  §5.1).  Also,  because  the  pdf  is  discrete,  we  calculate  the  mo¬ 
ments  using  the  biased  estimator  Cy^y.  ~  ~J2r^r.iYr}j,  instead  of  the  usual 
unbiased  estimator  Cyty?  ~  ^-j-  ]T)r  Yr  ,Yrj.  The  initial  fields  and  pdf  are 
shown  in  Fig.  [2], 
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6.1.2.  Lock-Exchange  Verification  Benchmark:  Results  and  Discussion 

The  outputs  from  the  stochastic  run  are  reported  in  Fig.  [3]  for  both 
Grashof  numbers.  Comparisons  with  the  deterministic  runs  are  shown  in  Fig. 
[4]  and  Fig.  [5]  for  Gr  —  4  x  104  and  Gr  =  1.25  x  106,  respectively.  Finally, 
the  evolution  of  the  differences  between  the  stochastic  and  deterministic  runs 
for  both  Grashof  numbers  are  shown  in  Fig.  [6]. 

We  see  excellent  agreement  between  the  stochastic  realizations  and  the 
deterministic  runs  (Fig.  [5] -Fig.  [6])  for  this  challenging  benchmark.  Partic¬ 
ularly,  for  the  lower  Grashof  number  flow,  the  local  error  is  less  than  0.2% 
everywhere.  It  is  non-trivia!  that  the  complex  DO  implementation  is  capa¬ 
ble  of  reproducing  multiple  deterministic  runs  in  a  single  simulation.  Thus, 
based  on  these  results  and  many  other  tests  (not  shown),  we  conclude  that 
our  implementation  is  correct,  that  is,  we  are  solving  the  intended  equations. 

The  growth  of  differences  between  the  deterministic  and  stochastic  sim¬ 
ulations  (Fig.  [6])  should  be  explained  by  the  main  differences  between  the 
stochastic  and  deterministic  solvers,  in  particular  the  advection  schemes,  and 
the  evolution  of  the  stochastic  coefficients.  We  found  that  the  magnitude  of 
the  differences  over  time  is  larger  for  coarser  space  and  time  resolution  runs. 
This  suggests  the  error  is  due  to  spatial  and/or  temporal  truncation  error. 
The  advection  scheme  does  not  contribute  significantly  to  the  error  at  the  re¬ 
ported  resolution,  since  using  a  CDS  advection  scheme  instead  of  the  TVD* 
scheme  for  the  stochastic  modes  did  not  change  the  reported  results  signifi¬ 
cantly.  Reducing  the  time-step  to  At  =  1/1024  reduced  the  error  at  the  final 
time  from  ~2.1%  to  ~1.05%  for  the  higher  Grashof  number  flow.  This  indi¬ 
cates  that  the  error  is  dominated  by  a  temporal  truncation  compounded  error 
of  approximtely  0(At )  (as  expected,  see  §6.3).  The  primary  source  of  this 
compounded  error  may  be  from  the  evolution  of  the  stochastic  coefficients 
and/or  the  modes. 

We  examine  the  errors  more  closely  to  determine  the  primary  source  of 
the  small  numerical  errors.  Note  that  the  density  differences  are  either  a  bit 
too  small  (first  realization  in  Fig.  [5])  or  a  bit  too  large  (last  three  realizations 
in  Fig.  [5]),  causing  phase  errors  as  observed  around  the  density  interface. 
These  phase  errors  can  be  caused  by  small  errors  in  the  relative  magnitudes  of 
the  stochastic  coefficients.  We  found  that  different  orthonormalization  strate¬ 
gies  (§5.2),  which  affect  the  relative  magnitudes  of  the  stochastic  coefficients, 
can  have  also  impact  the  magnitude  of  the  phase  errors.  In  particular,  when 
using  a  Gram-Schmidt  orthonormalization,  the  average  difference  was  greater 
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Mean  Mode  1  Mode  2 


Streamlines 


Figure  3:  The  DO  mean,  modes  at  non-dimensional  time  t  =  5,  and  evolution  of  stochastic 
coefficients.  For  Gr  =  1.25  x  106  (top)  and  Gr  =  4x  104  (bottom).  The  higher-Gr  flow 
has  sharper  gradients,  and  its  coefficients  are  larger  than  the  lower-Gr  flow.  Streamlines 
shown  over  density  in  color. 
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Final  Time:  Gr  =  1 ,25e6  Final  Time:  Gr  =  4e4 


Normalized 

difference 


Figure  4:  The  four  DO  realizations  (top  row),  the  deterministic  runs  (middle  row),  and 
the  DO  realizations  minus  the  deterministic  runs  normalized  by  ||  ^deterministic  ||  2  (bottom 
row)  for  Gr  =  4  x  104  (resolution  256x256  with  512  x  5  time-steps).  The  stochastic 
solver  result  agrees  with  the  deterministic  solver  results,  which  validates  the  new  solution 
method.  Streamlines  shown  over  density  in  color. 
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Figure  5:  As  Fig.  [4]  but  with  Gr  =  1.25  x  106.  The  higher-Gr  flow  with  sharper  gradients 
has  a  larger  error  than  the  lower-Gr  flow. 


Figure  6:  The  error  for  both  Gr  flows  tend  to  grow  over  time  but  can  decrease.  In  both 
cases,  the  DO  mean  field  has  a  smaller  error  than  the  average  error  of  the  realizations. 
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than  3%  (or  50%  worse)  for  the  high  Gr  case  at  the  final  time.  Also,  consid¬ 
ering  that  the  larger  Gr  number  flow  has  larger  stochastic  coefficients  (Fig. 
[3])  than  the  lower  Gr  flow,  the  same  relative  error  in  the  stochastic  coeffi¬ 
cients  would  lead  to  a  larger  phase  error  (observed  when  comparing  Fig.  [4] 
to  Fig.  [5]).  Thus,  this  benchmark’s  results  indicate  that  the  primary  source 
of  error  originates  from  the  time  evolution  of  the  stochastic  coefficients. 

In  summary,  using  the  proposed  benchmark,  we  have  verified  the  correct 
implementation  of  the  stochastic  code.  Also,  based  on  the  results,  we  suggest 
that  the  error  is  dominated  by  the  temporal  truncation  error  of  the  evolution 
equations  for  the  stochastic  coefficients.  This  indicates  that  the  accuracy  of 
the  scheme  will  benefit  most  from  improving  the  temporal  discretization  for 
problems  that  require  long  time  integration. 

6.2.  Effect  of  Advection  Scheme 

The  purpose  of  this  benchmark  is  to  test  the  three  different  advection 
schemes  proposed  (§4.1.2).  We  again  use  a  version  of  the  lock-exchange 
problem  because  the  sharp  density  interface  will  highlight  numerical  oscil¬ 
lations.  We  modify  the  problem  by  introducing  symmetry,  which  should 
be  maintained  numerically.  Finally,  for  simplicity  we  only  consider  a  single 
stochastic  mode  with  a  bimodal  continuous  pdf. 

6.2.1.  Effect  of  Advection  Scheme:  Setup 

General  setup:  The  Schmidt  number  is  kept  constant,  Sc  =  1,  and  we 
present  results  for  Gr  =  1.25  x  106.  Initially  the  velocity  is  zero  everywhere, 
and  free-slip  boundary  conditions  are  used  at  the  boundaries  of  the  domain 
(Fig.  [1]).  The  domain  is  discretized  using  Ax  =  Ay  =  1/64,  and  At  = 
1/256.  A  lower  resolution  is  used  here  compared  to  the  cases  in  §6.1  in  order 
to  highlight  the  symmetry  errors  and  numerical  oscillations. 

Mean  initialization:  The  mean  density,  pressure,  and  velocity  are  all  zero 
everywhere,  initially. 

Mode  initialization:  The  density  profile  for  the  mode  is  the  same  normal¬ 
ized  hyperbolic  tan  profile  used  in  §6.1.  The  pressure  and  velocity  for  this 
mode  are  zero  everywhere,  initially. 

Stochastic  coefficient  initialization-.  The  bimodal  Gaussian  continuous  pdf 
is  represented  by  10,000  samples.  To  ensure  symmetry  in  the  problem,  the 
pdf  of  the  mode  must  be  carefully  specified.  We  initialized  the  pdf  by  first 
generating  2,500  samples  (W250o)  from  a  zero  mean  Gaussian  distribution  with 
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Figure  7:  Symmetric  lock-exchange  problem  ( Gr  =  1.25  x  106)  with  grid  resolution  64  x  64 
and  At  =  1/256  using  various  advection  schemes  for  the  modes.  Only  the  stochastic  den¬ 
sity  is  non-zero  initially,  with  a  bi-modal  pdf(top  row) .  The  two  most  extreme  realizations 
(largest  and  smallest  stochastic  coefficients)  are  plotted  in  each  case  (right  column).  Our 
new  averaged  TVD*  scheme  only  has  minor  oscillations  and  retains  symmetry  (third  row), 
while  the  CDS  advection  scheme  suffers  from  large  oscillations  (white  dashed  circles,  sec¬ 
ond  row),  and  the  one-sided  TVD  scheme  loses  symmetry  (red  dashed  circles,  last  row). 
Streamlines  shown  over  density  in  color. 
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standard  deviation  a  =  O.Ole1  ~  0.027.  Next  these  samples  are  duplicated 
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to  yield  10,000  samples  that  need  to  be  evolved.  These  initial  conditions  can 
be  seen  in  the  first  row  of  Fig.  [7]. 


6.2.2.  Effect  of  Advection  Scheme:  Results 

The  result  of  this  simulation  is  shown  in  Fig.  [7]  for  the  various  advection 
schemes.  While  the  CDS  scheme  maintains  symmetry  of  the  mean,  modes, 
pdf,  and  realizations,  some  clear  numerical  oscillations  are  present,  particu¬ 
larly  evident  in  the  realizations.  While  there  are  no  oscillations  in  the  TVD 
scheme,  it  clearly  loses  symmetry  in  the  mean,  modes,  pdf,  and  realizations, 
as  can  be  seen  (with  aid  of  the  dashed  guide  lines)  in  Fig.  [7].  The  TVD* 
scheme  only  develops  minor  oscillations,  which  can  be  barely  detected  when 
examining  the  realizations,  and  completely  retains  symmetry.  Thus,  the  new 
TVD*  scheme  is  the  preferred  scheme  among  the  four. 

Initially,  we  can  represent  all  density  realizations  exactly  with  only  one 
mode.  However,  as  different  densities  evolve  at  different  rates,  one  mode  is 
not  sufficient  to  represent  this  uncertainty.  Hence,  spurious  gradients  appear 
in  the  reconstructed  realizations.  We  purposely  chose  only  one  mode  to 
illustrate  this  limitation  of  methods  that  would  keep  the  number  of  modes 
fixed,  which  we  normally  do  not  [53]. 

In  summary,  we  have  concluded  that  the  TVD*  scheme  performs  ade¬ 
quately,  and  we  demonstrated  that  the  DO  implementation  can  reproduce 
vastly  different  realizations  for  a  problem. 


6.3.  Numerical  Convergence  Analysis 

The  purpose  of  this  benchmark  is  mainly  to  show  that  the  implemented 
scheme  is  converging.  Here  we  use  the  classical  lid-driven  cavity  flow,  and 
examine  the  numerical  convergence  under  spatial  and  temporal  refinement  of 
each  component  separately.  This  benchmark  does  not  have  a  variable  density, 
and  so  we  report  the  Reynolds  number  instead  of  the  Grashof  number.  We 
also  completed  convergence  tests  with  density,  with  results  analogous  to  those 
presented  below. 
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x  =  (0,1) 


u  =  (l,0R 


Figure  8:  The  lid-driven  cavity  flow  is  a  classical  benchmark  used  to  verify  convergence 
of  numerical  implementations.  Uncertainty  for  this  case  will  be  introduced  through  the 
initial  conditions. 

6.3.1.  Numerical  Convergence  Analysis:  Setup 

General  setup:  We  present  results  for  a  Reynolds  number  of  Re  =  500, 
although  other  Reynolds  number  cases  (Re  G  [100,1000])  were  also  stud¬ 
ied,  giving  similar  results.  The  flow  is  driven  by  a  deterministic  bound¬ 
ary  condition  at  the  top  of  the  enclosed  box  Fig.  [8],  with  no-slip  velocity 
boundary  conditions,  and  uncertain  initial  conditions.  To  perform  the  con¬ 
vergence  study,  a  reference  solution  is  created  using  Ax  =  Ay  =  1/512,  and 
At  =  1/4096,  which  is  sufficient  resolution  at  this  Reynolds  number  [12,  11], 
We  examine  the  difference  between  the  reference  solution,  and  simulations 
with  coarser  space  and  time  resolutions  by  interpolating  (using  splines)  the 
fine  solution  onto  the  coarse  resolution  grid,  and  taking  the  L 2  norm  over  the 
interior  of  the  domain,  T>j  G  [0.25,0.75]  x  [0.25,0.75],  to  avoid  the  boundary 
condition  singularities  at  the  top  two  corners. 

Mean  initialization:  For  a  challenging  case,  the  mean  velocity  and  pres¬ 
sure  are  initially  zero,  everywhere. 

Mode  initialization:  The  velocity  modes  are  initialized  by  specifying  the 
stream  function 

^m,n(x,  y )  =  Cm,n  sin(7nc)  sin(7rMx)  sin(7n/)  sin(nNy), 

where  Cm,n  =  y  \J N (f ^  (f )5^M  ^  is  the  normalization  con¬ 

stant;  the  delta  function  5(x )  takes  the  value  1  if  x  =  0  and  0  otherwise.  The 
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velocity  modes  are  then  specified  as 


Ui  = 


d_ 

dy 


V>( M,N)t 7 


Vi  ~  dx^(M’N)i: 


where  (. M,N )  =  {(1, 1),  (1,2),  (1,3)}.  The  initial  pressure  for  the  modes  is 
specified  as  zero  everywhere. 

Stochastic  coefficient  initialization:  The  pdf  is  created  using  5,000  samples 
of  zero  mean  Gaussian  distributions  with  variances  Var(Kr5000jj)  = 

The  Yt  samples  are  then  created  using  a  similar  procedure  used  in  §6.2 


Y*  = 

1  r,i 


Yr 


r  5000 


~Yr 


r5000 


To  ensure  that  the  final  generated  samples  have  numerical  variances  exactly 
as  specified,  we  correct  the  samples  based  on  the  numerically  calculated  vari¬ 
ance 


t  -  _ w*  y/ Vsr  ( 1^5000  ,j )' 

r,i  *  ’ 

where  Var9(arjj)  =  yA-  Yll=i  ar,i  is  Hie  calculated  sample  variance.  The  ini¬ 
tialization  for  this  problem  can  be  seen  in  the  first  row  of  Fig.  [9]. 


6.3.2.  Numerical  Convergence  Analysis:  Results  and  Discussion 

Three  time  snapshots  of  the  reference  solution  are  shown  in  Fig.  [9]. 
While  the  marginal  seems  to  remain  approximately  Gaussian,  the  sample- 
scatter  diagram  clearly  shows  the  very  non-Gaussian  behavior  for  this  bench¬ 
mark.  Also,  we  confirm  from  Fig.  [10]  that  the  variances  of  the  modes  are 
decreasing,  which  is  expected  since  this  problem  has  a  deterministic  steady- 
state  solution.  Thus,  the  stochastic  solution  is  behaving  as  expected. 

Next,  examining  the  convergence,  we  see  that  the  numerical  error  for 
all  components  decrease  with  temporal  and  spatial  refinement.  The  conver¬ 
gence  is  near  optimal  at  large  grid  sizes  for  all  variables.  These  results  with 
the  velocity  components  separated  are  tabulated  in  Table  [1]  and  Table  [2], 
which  also  shows  that  the  stochastic  coefficients  are  converging  optimally. 
Even  though  a  fourth-order  RK  method  is  used  to  advance  the  stochastic 
coefficients,  the  total-DO  convergence  is  first  order  (based  on  choices  made 
in  §3.5).  Thus,  we  observe  near-optimal  convergence  for  all  variables,  which 
suggests  that  the  implementation  is  correct. 
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Figure  9:  Evolution  of  reference  solution  for  the  lid-driven  cavity  flow  ( Re  =  500)  with 
stochastic  initial  conditions  over  time  (resolution  512x512  with  4096  x  5  time-steps). 
This  case  demonstrates  that  the  DO  representation  with  our  implementation  is  able  to 
capture  and  evolve  non-Gaussian  statistics  (see  sample  scatter  plot)  for  a  continuous  pdf. 
The  marginal  pdfs  are  normalized  on  the  plot,  and  the  variance  can  be  read  off  Fig.  [10]. 
Streamlines  shown  over  pressure  in  color.  The  sample  scatter  is  colored  by  the  s  =  3 
stochastic  coefficient. 
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Figure  10:  Evolution  of  the  stochastic  energy  (Var(Y)))  for  the  lid-driven  cavity  reference 
solution.  The  tick-marks  on  the  time-axis  corresponds  to  the  time  snapshots  in  Fig.  [9]. 


At  Ax 


Figure  11:  The  control  volume  size  is  held  fixed  at  Ax  =  1/512  for  the  time  convergence 
(left),  and  the  time  step  size  is  held  fixed  at  At  =  1/4096  for  the  spatial  convergence 
(right).  The  error  ( 1 1 </>ref  —  4>n || 2)  decreases  with  both  temporal-  and  spatial-refinement 
for  each  component,  and  convergence  is  near-optimal  (order  1  in  time  and  2  in  space). 
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L6 

1.2 


Mode  2 
||  ^  ||  2 
l.le-03 
3.2e-03 
5.7e-04 
1.7e-03 
6.4e-04 
1.9e-03 
1.8e-04 
5.4e-04 


O 

1.6 

1.2 

L6 

1.2 

L6 

1.2 

L6 

1.2 


Mode  3 
||  ^  ||  2 
8.4e-04 
2.5e-03 
7.8e-04 
2.3e-03 
8.6e-04 
2.5e-03 
3.6e-04 
l.le-03 


O 

1.6 

1.2 

L6 

1.2 

L6 

1.2 

L6 

1.2 


Table  1:  Temporal  convergence  of  lid-driven  cavity  flow.  Tabulated  is  the  error  (e  = 
1 1 </>ref  —  0jv*  II2)  between  the  reference  solution  and  the  solution  using  At  =  1/Nt,  and  the 
approximate  order  of  convergence  O.  This  is  tabulated  for  the  grid  size,  Ax  =  1/512. 


6.4.  Stochastic  Convergence 

The  purpose  of  this  benchmark  is  to  further  assess  the  performance  of  the 
code,  demonstrate  that  a  large  number  of  modes  can  be  used,  and  study  the 
effect  of  the  stochastic  discretization.  Here,  we  extend  the  classic  shedding  of 
vortices  by  a  uniform  flow  as  it  encounters  a  symmetric  obstacle  to  stochastic 
DO  computations.  Since  this  problem  is  symmetric,  the  stochastic  solution 
should  also  be  symmetric.  Numerical  or  artificially  introduced  perturbations 
initiate  the  non-symmetric  laminar  shedding  of  vortices.  However,  if  these 
perturbations  are  symmetric,  there  should  be  no  preferential  direction  for 
vortex  shedding.  Thus,  we  expect  that  a  carefully  initialized  simulation  of 
the  stochastic  solver  should  be  able  to  capture  both  shedding  directions,  and 
this  expectation  presents  an  excellent  opportunity  to  assess  performance.  We 
also  quantify  the  effects  on  accuracy  of  the  number  of  stochastic  samples,  and 
their  time-order  of  integration. 

6.4.I.  Stochastic  Convergence:  Setup 

The  benchmark  is  an  open  flow  in  a  frictionless  pipe  with  a  square  cylin¬ 
drical  obstacle  (Fig.  [12]),  which  is  a  classic  test  for  deterministic  flow  solvers. 

General  setup:  We  present  results  for  a  Reynolds  number  of  Re  =  100, 
although  other  cases  were  also  studied.  The  flow  is  driven  by  a  determin¬ 
istic  uniform  inlet  boundary  condition  (left  of  domain),  with  slip  velocity 
boundary  conditions  at  the  top  and  bottom,  open  boundary  conditions  at 
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•  Nx 

Mean 
||  6 1|  2 

e> 

Mode  1 

Iklls 

£> 

Mode  2 
||  ^  ||  2 

£> 

Mode  3 
||  6 1|  2 

£> 

P 

256 

4.0e-05 

2.5 

1.5e-04 

2.2 

3.7e-04 

2.6 

6.9e-04 

2.3 

128 

2.2e-04 

2.5 

6.9e-04 

2.0 

2.3e-03 

2.4 

3.5e-03 

2.0 

64 

1.2e-03 

2.3 

2.7e-03 

2.1 

1.2e-02 

2.3 

1.4e-02 

1.5 

u 

256 

4.7e-05 

2.4 

1.2e-04 

2.4 

3.1e-04 

2.2 

8.3e-04 

2.4 

128 

2.4e-04 

2.4 

6.4e-04 

2.4 

1.5e-03 

2.1 

4.4e-03 

2.3 

64 

1.3e-03 

2.2 

3.3e-03 

2.3 

6.2e-03 

2.2 

2.1e-02 

2.0 

V 

256 

4.5e-05 

2.4 

1.3e-04 

2.4 

3.3e-04 

2.2 

9.2e-04 

2.4 

128 

2.4e-04 

2.2 

6.9e-04 

2.3 

1.6e-03 

2.1 

4.9e-03 

2.3 

64 

l.le-03 

2.2 

3.4e-03 

2.3 

6.5e-03 

2.2 

2.4e-02 

2.1 

Yi 

256 

l.le-04 

2.4 

1.6e-04 

2.4 

2.9e-04 

2.4 

128 

5.9e-04 

2.2 

8.7e-04 

2.3 

1.6e-03 

2.3 

64 

2.7e-03 

1.9 

4.2e-03 

2.1 

8.1e-03 

1.8 

Table  2:  Spatial  convergence  of  lid-driven  cavity  flow.  Tabulated  is  the  error  (e  =  ||dref  — 
<j)Nx  1 1 2 )  between  the  finest  (512  x  512)  and  present  (Nx  x  Nx)  grid,  and  the  approximate 
order  of  convergence  O.  This  is  tabulated  for  the  smallest  time  step,  At  =  1/4096. 


Figure  12:  Laminar  vortex  shedding  over  a  square  cylinder  in  a  channel.  Here  uncertainty 
originates  from  the  initial  conditions  and  uncertain  vortex  shedding:  depending  on  the 
perturbation,  the  first  vortex  could  either  be  shed  above  or  below  the  cylinder.  Within 
a  stochastic  framework,  however,  if  uncertainties  are  initially  symmetric,  the  mean  and 
modes  should  remain  symmetric,  and  this  is  evaluated  with  our  DO  numerics. 


38 


the  outlet,  and  symmetric  uncertain  initial  conditions.  All  simulations  use  a 
resolution  of  336  x  63  in  space,  and  63  x  40  in  time.  We  choose  to  integrate 
until  t  —  40  because  this  allows  the  statistics  to  reach  steady  values.  At 
t  =  40  the  mean  velocity  has  traveled  through  the  domain  2.5  times. 

Mean  initialization:  The  mean  velocity  and  pressure  are  initially  zero, 
everywhere. 

Mode  initialization:  While  the  exact  shape  of  the  initial  stochastic  per¬ 
turbations  are  not  important,  since  they  are  advected  out  of  the  domain,  it 
is  important  that  the  perturbations  are  symmetric.  We  initialize  the  velocity 
modes  by  specifying  the  stream  function 


fpM,N(x,y)  =  CM,Nsin{'Kx/a)sin{nMx/a)sin{/ny /b)sin(nNy /b)H 


where  Cm,n  is  the  normalization  constant  (as  in  §6.3),  and  a  =  16,  b  =  3 
are  the  width  and  height  of  the  domain  respectively.  The  velocity  modes  are 
then  specified  as 


d 

ui  = 


d 

Vi  = 


where 

(M,  N)  =  {(1, 1),  (2, 1),  (1,  2),  (3, 1),  (1,  3),  (2,  2),  (4, 1),  (1, 4),  (3, 2),  (2,  3)}, 

and  Bm  is  a  smoothing  function  created  numerically  from  the  domain  mask. 
Bm  is  created  from  the  mask  by  iteratively  averaging  each  control  volume  by 
its  own  value  and  its  four  neighbors  for  -A  iterations.  That  is,  at  iteration  k 

=  i  +  B^\i  -  l,j)  +  -  1) 

+h?M  1(*  +  1  ij)  +  -^M  1(hJ  +  1))  • 

The  initial  pressure  for  the  modes  is  specified  as  zero  everywhere. 

Stochastic  coefficient  initialization:  To  ensure  initial  symmetry,  the  sam¬ 
ples  for  the  stochastic  coefficients  are  created  using  the  same  procedure  as  in 
§6.3,  using  the  variances  V ar(Yrf  =  e(2-Mi-iVi)  (note  the  difference  of  +1  in 
the  exponential  from  §6.3).  The  reference  solution  uses  105  samples  for  the 
stochastic  coefficients,  and  a  4th  order  RK  time  integration  scheme  (§3.5). 
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Marginal  PDFs 


Figure  13:  Mean  field  and  first  5  modes  with  marginal  pdfs  at  final  non-dimensional  time 
T=40  for  Re=100,  and  the  evolution  of  the  mean,  (u,  u)^,  and  stochastic  energy,  Var(Y)), 
for  the  reference  solution  (resolution  63  x  336  with  63  x  40  time-steps) .  Streamlines  shown 
over  pressure  in  color..  Our  scheme  and  implementation  is  able  to  retain  symmetry  for 
the  most  important  first  four  modes. 
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Figure  14:  As  in  Fig.  [13],  but  showing  two  realizations  where  the  vortex  is  shed  in 
opposite  directions.  The  colorbar  for  the  pressure  is  as  on  Fig.  [13]. 

6-4-2.  Stochastic  Convergence:  Results  and  Discussion 

For  the  reference  solution,  we  observe  that  our  scheme  maintains  excel¬ 
lent  symmetry  for  the  most  important  first  four  stochastic  modes  (Fig.  [13]). 
From  the  realizations  (Fig.  [14]),  the  scheme  clearly  captures  both  shedding 
directions.  We  find  that  for  fewer  samples  and  lower  time  integration  ac¬ 
curacy,  the  symmetry  is  not  as  well  maintained  (not  shown).  This  can  be 
partially  explained  by  the  non-symmetric  bimodal  marginal  of  the  first  mode 
(Fig.  [16])  when  using  fewer  samples.  We  also  found  it  is  more  challenging 
to  maintain  symmetry  for  higher  Reynolds  numbers.  Nonetheless,  for  a  suffi¬ 
cient  number  of  samples  for  the  stochastic  coefficients,  our  scheme  maintains 
excellent  symmetry  for  all  Reynolds  numbers  we  tested. 

Also  from  Fig.  [13],  we  observe  that  the  variances  seem  to  reach  a  steady 
value  after  an  initial  transient  period.  The  smaller  variances  take  longer  to 
reach  a  steady  value,  with  the  highest  modes  perhaps  still  evolving  after  the 
final  time  step. 

For  our  implementation,  the  reference  DO  simulation  (Fig.  [13])  was  com¬ 
puted  in  approximately  3  hours  a  2.4Ghz  computer.  Each  time-step  required 
the  inversion  of  11  pressure  equations,  as  opposed  to  the  111  required  for  the 
scheme  without  the  stochastic  pseudo  pressure.  This  roughly  translates  to  a 
1000  %  increase  in  efficiency,  or  about  14  days  saving  in  terms  of  computa¬ 
tional  time  for  this  computer.  Since  then,  we  have  used  as  many  as  s  =  25 
modes  for  data  assimilation  applications  [58,  59]. 

Sample  sizes:  Examining  the  evolution  of  the  stochastic  coefficients  for 
different  sample  sizes  (top  of  Fig.  [15]),  we  see  that  the  differences  are  larger 
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Eneargy  Energy 


Figure  15:  Stochastic  energy  (Var(Y)))  of  the  stochastic  coefficients  over  time  for  different 
number  of  samples  (0( 4)  time  integration)  (top),  and  different  time  integration  schemes 
(10,000  samples)  (bottom).  The  overall  trends  are  well  captured  in  all  cases,  but  there 
are  noticeable  errors  for  less  energetic  modes  after  long  integration  times  when  a  small 
number  of  samples  and  a  low  order  time  discretization  scheme  is  used.  A  relatively  small 
number  of  samples  (10,000)  can  be  used  for  this  benchmark  since  nearly  identical  results 
compared  to  100,000  were  found  for  the  most  energetic  four  modes. 
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Samples:  1e2 
Order:  0(4) 


Mode  1 


1e3 

0(4) 


Mode  2 


Mode  3 


Mode  4 


Mode  5 


1e4 


0(1) 


1e5 

0(4) 


Figure  16:  Marginal  probability  density  functions  of  first  five  modes  for  the  flow  over 
a  square  cylinder  benchmark  at  the  final  time  (columns  1-3,  increasing  sample  sizes; 
columns  4-5,  decreasing  order  in  time;  column  6,  reference).  Using  10,000  samples  and 
high  order  integration  for  the  Stochastic  coefficients,  the  continuous  marginal  pdfs  are  well- 
represented  by  our  approach,  although  the  bimodal  peaks  are  not  symmetric.  Compared 
to  the  fourth  order  time  integration  scheme,  the  second  order  time  integration  scheme 
has  similar  marginals  for  the  first  5  modes,  whereas  the  first  order  scheme  differs  after 
the  second  mode.  The  marginals  have  similar  shapes  for  all  sample  sizes,  but  the  best 
representations  have  larger  sample  sizes. 
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for  higher  modes.  For  the  most  energetic  first  four  coefficients,  the  results 
agree  well  with  the  reference  solution  when  using  >  1,000  samples.  After 
the  5th  coefficient,  the  differences  become  larger,  and  in  the  10th  coefficient, 
the  differences  begin  sooner  and  have  larger  relative  magnitudes.  Note  that 
the  logarithmic  scale  magnifies  smaller  errors,  and  the  differences  for  the 
10th  are  fully  insignificant  when  compared  to  the  variance  of  the  first  mode. 
Thus,  with  a  relatively  small  number  of  samples,  it  is  possible  to  capture  the 
stochastic  coefficient’s  variances  accurately.  Also,  it  appears  as  though  the 
solution  converges  when  the  number  of  samples  increases,  which  increases 
our  confidence  in  this  approach. 

Not  only  are  the  evolution  of  the  variances  well-represented  by  a  smaller 
number  of  samples,  the  shape  of  the  marginals  are  also  well- reproduced  (Fig. 
[16]).  We  observe  that  the  marginals  are  well-reproduced  using  >  1,000 
samples,  and  for  all  time  integration  schemes  (Fig.  [16]).  The  magnitudes 
of  the  bimodal  peaks  in  the  first  mode  are  more  symmetric  with  increased 
sample  sizes,  which  suggests  using  upwards  of  10,000  samples  is,  perhaps, 
advised.  In  fact,  using  such  a  number  of  samples  does  not  affect  the  overall 
cost  of  the  solution  scheme.  That  the  marginals  are  well-represented  using 
smaller  sizes  is  encouraging,  since  it  suggests  that  a  small  number  of  samples 
(that  won’t  negatively  impact  the  efficiency  of  the  scheme)  could  be  used  for 
problems  with  large  s:  this  would  not  substantially  affect  accuracy. 

Order  in  time:  Recall,  we  are  only  changing  the  order  of  the  ODE  solver 
used  to  evolve  the  stochastic  coefficients,  but  the  overall  order  of  the  scheme  is 
kept  fixed  at  first  order  (see  §3.4).  Examining  the  evolution  of  the  stochastic 
coefficients  for  different  time  integration  schemes  (bottom  of  Fig.  [15]),  we 
see  that  the  second  and  fourth  order  schemes  agree  well  for  all  coefficients. 
These  results  indicate  that  a  first  order  time  integration  scheme  can  results 
in  significant  errors,  and  should  be  avoided. 

In  summary,  we  have  simulated  a  stochastic  version  of  the  classical  flow 
over  a  square  cylinder  benchmark.  We  found  that  our  stochastic  solution 
maintains  excellent  symmetry,  and  our  new  scheme  is  1000%  more  efficient 
when  using  10  stochastic  modes.  We  also  found  that  the  statistics  seem  to 
converge  with  increased  number  of  samples  and  that  accuracy  benefits  from 
using  a  second-order  ODE  solver,  even  though  the  overall  scheme  is  limited 
to  first  order.  Finally,  smaller  sample  sizes  of  the  stochastic  coefficients  still 
produced  accurate  time-evolving  statistic  and  marginal  pdfs,  which  is  very 
encouraging. 
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7.  Conclusions  and  Future  Research 

We  have  derived  efficient  computational  schemes  for  the  DO  methodol¬ 
ogy  applied  to  unsteady  stochastic  Navier-Stokes  and  Boussinesq  equations, 
and  illustrated  and  studied  the  numerical  aspects  of  these  schemes.  For 
the  discretizations  in  time,  we  developed  semi-implicit  projection  methods 
for  the  mean  and  the  modes,  and  we  employed  time-marching  schemes  of 
first  to  fourth  order  for  the  stochastic  coefficients.  For  the  discretizations  of 
the  physical  space,  we  employ  conservative  second-order  finite-volumes,  with 
a  special  treatment  for  the  advection  terms  based  on  a  TVD  scheme  with 
monotonized  central  symmetric  flux  limiter.  Several  alternate  discretizations 
in  time  and  space  were  outlined  and  illustrated  by  comparisons.  We  also  ad¬ 
dressed  several  numerical  issues  specific  to  the  DO  method  for  fluid  and  ocean 
flows.  In  particular  we  have  shown:  how  to  define  pseudo-stochastic  pres¬ 
sures  to  reduce  the  number  of  matrix  inversions  for  the  pressure  from  0(s2) 
to  O(s)]  how  to  treat  advection  by  the  stochastic  modes  using  symmetric 
approximate  TVD  schemes;  how  to  deal  with  singular  subspace  covariances 
by  generalized  inversion;  and  how  to  maintain  orthonormal  modes  at  the  nu¬ 
merical  level  and  so  account  for  truncation  and  round-off  errors  during  time 
integration.  Finally,  we  evaluated  our  schemes  using  a  varied  set  of  stochastic 
flows,  hence  illustrating  robustness  but  also  providing  benchmarks  for  future 
schemes  and  implementations. 

Using  an  asymmetric  Dirac-stochastic  lock-exchange  benchmark,  we  found 
excellent  agreement  among  multiple  realizations  from  a  deterministic  code 
and  the  realizations  generated  from  a  single  stochastic  simulation.  This 
validated  our  numerical  implementation,  and  confirmed  that  the  pseudo¬ 
stochastic  pressure  approach  works  in  practice.  Using  a  symmetric  version 
of  the  lock-exchange  benchmark,  we  showed  that  the  symmetric  TVD  advec¬ 
tion  scheme  (TVD*)  works  well,  while  the  CDS  scheme  suffers  from  numerical 
oscillations,  and  the  non-symmetric  TVD  scheme  loses  symmetry. 

Using  a  lid-driven  cavity  flow  with  uncertain  initial  conditions,  we  showed 
that  each  component  converged  near-optimally  under  both  time  and  space 
refinement.  Additionally,  we  showed  that  even  when  a  higher-order  time  in¬ 
tegrator  is  used  for  the  stochastic  coefficients,  their  accuracy  is  still  limited  to 
approximately  first  order,  as  discussed  in  §3.4.  This  benchmark  also  demon¬ 
strated  an  expected  decay  in  the  variance  of  the  stochastic  coefficients,  which 
were  shown  to  be  very  non-Gaussian. 

Finally,  using  a  stochastic  flow  past  a  square  cylinder  in  a  confined  chan- 
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nel,  we  showed  that  the  stochastic  coefficients  converged  with  increased  sam¬ 
ples  sizes  and  orders  of  time  integration.  We  found  that  using  >  10,  000 
samples  and  a  second  order  time-integrator  yielded  adequate  performance. 
This  benchmark  also  demonstrated  that  our  discretized  DO  methodology 
successfully  captures  both  vortex  shedding  directions,  resulting  in  a  fully 
symmetric  mean  held.  Also,  we  noted  that  using  our  newly  defined  pseudo¬ 
stochastic  pressure,  we  were  1000%  more  efficient  computationally  when  us¬ 
ing  10  modes. 

Possible  future  studies  include  the  efficient  extension  of  our  present  frame¬ 
work  to  hows  with  uncertain  parameters  in  the  governing  equations  and 
with  stochastic  forcing  at  the  boundary  and  in  the  interior.  Also,  our  re¬ 
sults  suggest  that  simulations  with  long  time  integration  will  benefit  from 
more  accurate  time-integration  scheme.  Then,  while  our  proposed  symmet¬ 
ric  TVD-based  advection  and  orthonormalization  schemes  for  the  stochastic 
modes  were  shown  to  perform  adequately,  more  optimal  treatments  are  pos¬ 
sible.  Alternative  means  of  discretizing  the  stochastic  coefficients  are  worth 
exploring  for  applications  where  rare  events  are  important.  Finally,  specihc 
multi-resolution  DO  schemes  can  be  derived  for  multiscale  stochastic  fluid 
and  ocean  hows.  Unstructured  grids  as  well  as  nested  approaches  [8,  20] 
will  then  be  useful.  The  utilization  of  the  present  schemes  as  well  as  these 
future  advances  will  allow  the  prediction  of  uncertainty  in  realistic  simula¬ 
tions  of  multi-physics  fluid  systems,  over  a  wide  range  of  applications,  from 
micro-nano  fluid  engineering  to  multiscale  ocean  and  climate  studies. 
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Appendix  A.  DO  Discretization 

In  this  Appendix  we  present  the  stochastic  DO  equations  for  Boussinesq 
dynamics,  and  briefly  summarize  the  required  steps  (from  [52,  51,  54])  to 
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obtain  these  equations. 

Starting  with  a  generalized  Karhunen-Loeve  expansion  [34,  51],  one  de¬ 
composes  the  solution  of  eqns.  (l)-(3)  into  a  Dynamically  Orthogonal  (DO) 
field  expansion  [52,  51,  53]  for  the  velocity,  density  and  pressure3  fields 


i= i 

$  =  $  + 

s  s  s 

p(x,  f;  u)  =  p(x,  f)  +  y,  Yi(t,  uj)pi(x,  t)  +  EE 

2=1  2=1  J  =  1 

V  =  P  +  E/h  +  Y^j'Piji 

(A.2) 


where  $  =  [$1,  <f>2, . .  ,]7  =  [u,  p]T  is  the  vector  of  prognostic  state  variables. 
The  scalar  s  =  s{t)  defines  the  time-dependent  dimension  of  the  stochas¬ 
tic  subspace,  i.e.  s  is  a  discrete  number  of  stochastic  terms  retained  from 
a  complete  expansion.  The  field  functions  3>,;(x,  t)  are  the  s  orthonormal 
deterministic  modes  and  the  Yi{t\uj)  are  their  s  zero-mean  stochastic  coeffi¬ 
cients,  in  general  non-Gaussian.  In  our  notation,  we  use  the  Einstein  sum¬ 
mation  exclusively  for  summations  related  to  the  stochastic  expansion.  The 
decomposition  of  pressure  into  a  mean,  linear  modal  and  quadratic  modal 
component  follows  from  the  Pressure  Poisson  Equations  (see  §3.2  or  Sapsis 
and  Lermusiaux  [52]).  Since  both  the  modes  and  stochastic  coefficients  are 
functions  of  time,  a  redundancy  arises,  which  is  resolved  by  the  DO  condition, 

=°,  Vi,  j  G  [1,  2, . . . ,  s],  (A. 3) 

where  the  inner  product  is  defined  as  (a,  b)^  =  jv  JA(a*&*)  cW  for  arbitrary 
vectors  of  spatial  functions  a  =  [a1,  a2, . .  .]T  and  b  =  [61,  b2, . .  ,]T. 

Using  the  DO  condition,  an  exact  set  of  equations  can  be  obtained  that 
governs  the  evolution  of  the  mean,  modes,  and  stochastic  coefficients  of  the 
generalized  Karhunen-Loeve  expansion.  The  only  approximation  arises  from 


3For  convenience,  we  changed  the  sign  of  pij  from  the  original  definition  in  the  Sapsis 
and  Lermusiaux  [52] 


47 


the  truncation  of  the  DO  expansion  to  s(t)  terms.  First  substitute  (A.l)  into 
eqns.  (l)-(3)  to  obtain  (using  a  Langevin  notation): 


wu  ain  u  u  j  , 

-Qj;  +  -^-Uj  +  F)—  =  C  (u  +  YiUi,p,x,t;uj) , 
dp  dYi  „ dpi 

m  +Hrp‘  +  Y'W  =  c 


(A. 4) 


and  DO  decomposed  versions  of  (2)- (3).  It  is  from  these  equations,  within 
which  the  DO  decomposition  was  inserted,  that  the  equations  for  the  mean, 
modes  and  their  coefficients  are  obtained,  using  the  expectation  operator, 
the  spatial  inner  product,  and  eqns.  (A.2)-(A.3). 

Mean.  To  obtain  a  rate  of  change  for  the  mean  fields,  the  idea  is  to  elim¬ 
inate  the  random  components  in  the  left-hand-sides  of  eqns.  (A. 4).  Hence, 
taking  the  expectation  of  eqns.  (A. 4)  it  can  be  found  that  the  evolution  of 
the  mean  fields  are  governed  by 


V  •  u  =  0,  x  G  V, 

W  ~  =  "v ' ~vp  +  peS 

-  C ViYj  ( Vpij  +  V  •  (uj-Ui)) ,  x  g  D. 
^ - ^V2P  =  -V  •  (up)  -  Cyy  V  •  (u.'Pi),  X6D, 

dt  ScVG^  3 

with  the  initial  and  boundary  conditions  given  by 


u  (x,  0)  =  u0,  x  G  V, 
p(x,0)  =  p0,  x6D, 


U  =  g  D, 

<9u 

=  Sjv, 
on 

P  =  9  Dpi 
dp 

m = 3N- 


x  g  dVD, 
x  G  dT>  jv, 
X  G  dVDp, 
x  G  dVNp, 


(A.5) 


(A.6) 


(A.  7) 


where  Cy4y.  =  Eu  \YiYj ]  is  an  element  of  the  covariance  matrix  in  the  stochas¬ 
tic/error  subspace.  Deterministic  initial  and  boundary  conditions  (i  quan¬ 
tities)  are  assigned  to  the  mean.  Note  that  in  general  the  vector  V  •  (UjUj) 
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differs  from  V  •  (ujUj)  (recall  that  V  •  (u,-Uj)9  =  <)u^Uiq  ^  ~W^JL e-g-, 

dvjUi  /  dviiij  \ 
dy  '  dy  '  * 

Modes.  The  evolution  of  the  modes  is  also  obtained  from  eqns.  (A. 4). 
To  do  so,  the  idea  is  to  eliminate  the  random  coefficients  in  front  of  the 
time  derivatives  of  the  modes.  The  essential  steps  are  to  multiply  these 
equations  with  a  stochastic  coefficient  Y),  apply  the  expectation  operator,  and 
substitute  an  expression  for  E ^  [^Y)] ,  which  is  obtained  by  projecting  the 
equation  unto  and  imposing  the  DO  condition.  From  this  the  following 
governing  evolution  equations  for  the  modes  can  be  found: 


V  •  Uj  =  0,  xtD, 

<9u 

=  Qj  —  (Qi,  *&j)v  ui  ,  x  ^  ^1  (A. 8) 

^  =  Qi  -  (Qi,  ®j)v  PE  x  e  v, 


where 


Qi 

Qr 


Qi 


[Qr,  Q?f  = 

u 

J  [£"  yt] 

,  Cy]Y.E“  [CpYj\ 

^  - 

V  •  (ujU 

)  -  V .  ( 

uuj)  -  Vpi  +  pp 

-  Cyi1yjMr.ymy„  (Vprun  +  V  • 

(unum) )  , 

ScVgE  p' 

> 

1 

> 

1 

la 

•  ip-pi) 

-  CAMyymy„  (V  • 

iy^nPm)  ) 

1 

and  Myjymyn  =  Eu  [YjYmYn]  is  a  third  moment.  The  right-hand-sides  in 
eqns.  (A. 8)  correspond  the  total  rate  of  change  of  the  subspace  (without  a 
DO  condition)  minus  the  projection  of  this  rate  of  change  on  the  subspace 
itself  (which  is  subtracted  to  ensure  the  DO  condition).  We  note  that  in 
general  V  •  (uUj)  ^  V  •  (u,u)  since  ^  for  example. 

The  initial  and  boundary  conditions  for  the  modes  are  obtained  from 
those  of  the  full  stochastic  fields,  eqns.  (2)  and  (3),  but  reduced  to  their 
dominant  initial  error  subspace  of  size  .s'0, 


u i  (x,  0)  =  Uj  0  (x) ,  x  G  V, 
Pi  (x,  0)  =  pifi  (x) ,  x  G  V, 


(A.9) 
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and 


<9u£ 

dn 

Pi 

dpi 

dn 


Si,D  ■ 

X  G  dVD, 

g  i,N, 

x  g  dv N% 

Pi, Dp  i 

X  G  dVDp 

Pi,Np, 

x  G  dVNp 

(KAO) 


Coefficients.  Finally,  to  obtain  the  evolution  of  the  stochastic  coefficients 
from  eqns.  (A. 4),  the  idea  is  to  eliminate  the  modes  in  the  term  contain¬ 
ing  the  time  derivatives  of  the  random  coefficients.  To  do  so,  project  the 
evolution  eqns.  (A. 4)  onto  each  mode  i,  apply  the  DO  conditions,  and  essen¬ 
tially  impose  that  each  coefficient  is  of  zero  mean.  The  resulting  governing 
stochastic  ODEs  are 


<m 

dt 


(Fm-  #i>o  Ym  -  {Vpm  +  V  ■  (UnUm),  Ui)v  (YmYn  -  Cy„yJ  (A-U) 
(V  •  (ll nPm)i  Pi)x>  mYn  C YmYn )  j 


where  £  =  [£u,  jCp]t  and  Fm  =  [F“  ,  Fpn]T  with 


Fu  =  — 

m  VG~r 


V2um,  —  V  ■  (umu)  —  V  •  fuur 


Vpm  +  pme9, 


F^ 

m 


V  •  (u mp)  ~  V  •  (u Pm); 


The  initial  conditions  for  the  coefficients  are  obtained  from  those  of  the 
full  stochastic  fields,  eqn.  (2),  by  projection  onto  each  initial  mode  i,  ^o, 
and  removal  of  the  mean.  This  leads  to: 


(t0;  u)  —  ($o  —  <&o  ,  j 

=  (U0  —  U0  ,  +  (p0  ~  Po  ,  Pi, o)V  )  (aJ  E  Fl  . 


(A.12) 


In  this  Appendix  we  stated  the  DO  equations  for  the  mean,  modes,  and 
stochastic  coefficients  of  the  stochastic  Boussinesq  equations. 
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a) 

Orthonormalize 
basis  without 
modifying 
coefficients 


IT=y|  v 

b> 

Orthonormalize 
basis  and  retain 
"O^1  ,Y2 )  realization 

energy  and 
^  'direction' 


Figure  B.17:  There  are  different  ways  to  treat  the  stochastic  coefficients  when  correcting 
the  discrete  orthonormalization  of  the  modes.  In  this  pictorial  analogy,  the  standard  2D 
Cartesian  bases  and  coordinates  represent  the  modes  and  stochastic  coefficients,  respec¬ 
tively. 


Appendix  B.  Numerical  orthonormalization  procedure 

In  this  Appendix,  we  present  our  newly  developed,  numerically  efficient, 
orthonormalization  procedure  based  on  the  discussion  in  §5.2. 

Consider  the  pictorial  analogy  showing  different  strategies  when  orthonor- 
malizing  the  basis  (Fig.  [B.17]).  In  case  (a),  the  argument  is  that  the  coef¬ 
ficients  were  evolved  correctly,  and  all  the  error  lies  in  the  basis  evolution. 
This  is  not  true  because  we  evolve  the  coefficients  by  taking  an  inner  prod¬ 
uct  of  the  equations  using  the  bases  at  the  old  time-step.  In  case  (b),  the 
argument  is  that  the  energy  and  ‘direction’  of  the  realizations  were  correctly 
evolved.  When  correcting  the  bases,  the  direction  and  energy  are  maintained, 
but  the  actual  realization  is  modified.  In  case  (c),  the  argument  is  that  the 
realizations  were  evolved  correctly.  However,  it  is  unlikely  that  the  errors 
in  the  modes  and  coefficients  cancel  during  the  time  evolution.  Also,  while 
not  evident  from  the  pictorial  analogy,  the  energy  of  the  coefficients  change 
in  case  (c).  This  is  because  the  length  of  the  basis  changes  when  orthonor- 
malized,  and  this  length  change  is  then  accounted  for  in  the  coefficients  in 
order  to  maintain  the  coordinate  of  the  realization.  The  case  (b)  strategy  is 
a  compromise  between  cases  (a)  and  (c),  and  it  results  in  the  lowest  error  for 
the  lock-exchange  cases  in  §6.1.  The  errors  were  the  largest  for  case  (a),  and 
the  errors  for  case  (c)  were  marginally  larger  compared  to  case  (b). 

To  be  clear,  case  (b)  is  shown  to  perform  best  for  the  specific  cases  in 
§6.1,  and  these  results  are  not  proven  to  be  universally  applicable.  Also,  we 
made  no  attempt  to  provide  an  exhaustive  list  of  all  the  orthonormalization 
strategies.  For  example,  arguing  that  the  coefficients  still  live  in  the  bases 
of  the  old  time  step,  one  can  consider  projecting  the  coefficients  from  this 
old  basis  onto  the  new  basis  before  orthonormalization.  Therefore,  the  issue 
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of  the  optimal  orthonormalization  strategy  is  left  open  to  future  research. 
Nonetheless,  we  believe  our  current  approach  is  sufficient  for  most  applica¬ 
tions  since  no  strategy  will  influence  the  scheme’s  overall  order  of  accuracy. 
Any  improved  strategy  will  only  reduce  the  size  of  the  error,  which  will  be  of 
order  At  smaller  than  the  truncation  error.  Therefore,  any  accuracy  issues 
clue  to  orthonormalization  can  be  resolved  by  reducing  the  size  of  the  time 
step  when  using  the  current  scheme.  Our  orthonormalization  procedure  for 
case  (b)  is  described  next. 

For  case  (b),  we  want  the  new  orthonormalized  basis  to:  approximately 
recover  the  specific  realizations  before  the  orthonormalization;  and  exactly 
maintain  the  stochastic  energy.  Mathematically  we,  ideally,  want  both: 


rO\T 


Tr(C 


k* 


=  Tr(C 


k*,0 


(B.l) 

(B.2) 


where  the  orthonormalized  basis  is  defined  such  that  (&h,n  ^ hj )v  —  and 

the  trace  operator  Tr  simply  sums  the  diagonal  entries  of  a  matrix.  In  order 
to  preserve  the  second  property,  we  first  rotate  the  stochastic  system  such 
that  the  covariance  matrix  is  diagonal: 


C‘y.ri  =  E[V, =  VdcDDcV50, 

Yjf  =  Y,.\'nc: 

*£?  =  *wVdc, 


(B.3) 

(B.4) 

(B.5) 


which  gives  the  decorrelated  (or  DC)  stochastic  coefficients.  By  starting 
with  a  decorrelated  system,  we  have  a  unique  and  convenient  reference  frame 
which  we  can  use  to  enforce  the  Tr(Cy?y )  property. 

Next  we  perform  the  orthonormalization  as  follows 


Mij 

yOC 


.  <■>£?}„  =  VmDmVI,, 

(B.6) 

Yrf\M^D^h 

(B.7) 

(B.8) 

The  above  could  be  accomplished  by  directly  calculating  the  Singular  Value 
Decomposition  of  the  matrix,  but  the  above  procedure  has  the  advan¬ 
tage  of  enabling  the  user  to  specify  any  desired  (possibly  non-linear)  inner 
product.  In  either  case,  the  two  approaches  are  similar  in  efficiency  and  give 
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the  same  results.  To  verify  that  (B.l)  is  satisfied,  we  substitute  the  above 
expression  in  (B.l) 


*£F  {yrf)T  = 


OC /\sOC\T 


Vdm  nfvMy D 


=  Vm  y  D^1  v/D^V  m 

=  $>?c(Ydc)t 

h,i  V  r,i  )  i 


M 
DC\T 


since  the  eigen- vectors  are  orthonormal  V mY\4  =  I.  Therefore  the  first 
property,  (B.l),  is  exactly  satisfied  at  this  stage.  Stopping  here  would  give 
the  case  (c)  strategy. 

Finally,  we  decorrelate  the  samples  again,  while  ensuring  that  the  stochas¬ 
tic  energy  is  preserved: 


C“°  =  E[  yf  yf]  =  VoDoVj, 


-OCVOC 


Y°  =  Y°CY , 


■  r,i 


r,i 


o\ 


'Tr(D 


DC) 


Tr(D 


o) 


K  = 


(B.9) 

(B.10) 

(B.ll) 


Note  that  Tr(Do)  =  ^(DdcDm),  but  that  Do  ^  D^cDm-  This  ensures 
that  (B.2)  is  maintained,  but  now  we  can  only  approximately  satisfy  (B.l). 
The  over-all  correction  is  then  as  follows: 


Y%  =  Yr,i Vdc  Vm \/D mVo  J  Tr^of  ’  (B'12) 

Ki  =  *h,iVDcVMyfDT}Vo.  (B.13) 

This  procedure  is  robust  even  if  C  1pYc  is  singular,  however  it  does  require 
a  non-singular  Mg  which  will  normally  be  the  case  if  modes  are  properly 
initialized.  Also,  the  number  of  floating  point  operations  is  dominated  by 
the  calculation  of  Mg  if  the  number  of  spatial  points  exceeds  the  number 
of  realizations.  Nonetheless,  the  overall  cost  of  the  orthonormalization  is 
inexpensive  compared  to  the  inversion  of  the  pressure  Poisson  equations. 
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